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ABSTRACT 

The old open cluster M67 is an ideal testbed for current cluster evolution models 
because of its dynamically evolved structure and rich stellar populations that show 
clear signs of interaction between stellar, binary and cluster evolution. Here we present 
the first truly direct A^-body model for M67, evolved from zero age to 4 Gyr taking full 
account of cluster dynamics as well as stellar and binary evolution. Our preferred model 
starts with 36 000 stars (12 000 single stars and 12 000 binaries) and a total mass of 
^s^' \ nearly 19OOOM0, placed in a Galactic tidal field at 8.0 kpc from the Galactic Centre. 

■ Our choices for the initial conditions and for the primordial binary population are 

' explained in detail. At 4 Gyr, the age of M67, the total mass has reduced to 2 OOOMq 

\f~^ , as a result of mass loss and stellar escapes. The mass and half-mass radius of luminous 

C3 ' stars in the cluster are a good match to observations although the model is more 

centrally concentrated than observations indicate. The stellar mass and luminosity 
Q^. functions are significantly flattened by preferential escape of low-mass stars. We find 

I ' that M67 is dynamically old enough that information about the initial mass function 

O \ is lost, both from the current luminosity function and from the current mass fraction 

in white dwarfs. 

^ ' The model contains 20 blue stragglers at 4 Gyr which is slightly less than the 28 

observed in M67. Nine are in binaries. The blue stragglers were formed by a variety of 
means and we find formation paths for the whole variety observed in M67. Both the 
k> primordial binary population and the dynamical cluster environment play an essential 

, role in shaping the population. A substantial population of short-period primordial 

' binaries (with periods less than a few days) is needed to explain the observed number 

of blue stragglers in M67. The evolution and properties of two thirds of the blue 
stragglers, including all found in binaries, have been altered by cluster dynamics and 
nearly half would not have formed at all outside the cluster environment. On the other 
hand, the cluster environment is also instrumental in destroying potential BSs from 
the primordial binary population, so that the total number is in fact slightly smaller 
than what would be expected from evolving the same binary stars in isolation. 

Key words: stellar dynamics — methods: N-body simulations — stars: evolution — 
blue stragglers- binaries: close — open clusters and associations: general 



1 INTRODUCTION 

Star clusters have long been recognized as important tools 
for understanding many astrophysical processes. As such 
they are the focus of many observational programs, both 
from the ground (e.g. Kalirai et al. 2003; Kafka et al. 2004) 
and space (e.g. Grindlay et al. 2001; Piotto et al. 2002; 
Richer et al. 2004). Dynamical modelling had its beginnings 
over four decades ago (von Hoerner 1960; Aarseth 1966; van 
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Albada 1968) but it is only recently that the models have 
reached a state where genuine comparison with cluster ob- 
servations is possible. This is the result of software advances 
that have improved realism and speed as well as a huge in- 
crease in hardware performance. The realm of globular clus- 
ters is still out of reach of realistic direct A'-body models. So 
it is open clusters that garner initial attention in the attempt 
to confront cluster models with observations and vice- versa. 
In particular, old open clusters are of interest because these 
are dynamically well-evolved and offer the chance to observe 
the effects of interaction between stellar, binary and cluster 
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evolution. A good example is M67 which contains an abnor- 
mally large number of blue straggler (BS) stars for its age 
and size (Ahumada & Lapasset 1995) , as well as a good pro- 
portion of X-ray sources (van den Berg et al. 2004). Both 
are indicators that the cluster environment has affected the 
evolution of the stars. M67 is also close enough to be well 
studied and within range of current A'^-body modelling ca- 
pabilities^ . The idea is that by matching the characteristics 
of a model cluster and the observed cluster we can learn 
about the initial conditions, binary population and dynam- 
ical evolution of the cluster. In doing this we must bear in 
mind the need to marry theoretical and observational data 
reduction techniques (Portegies Zwart et al. 2004). 

The interest in M67 and its BSs stems not only from 
the large number of BSs but also from their diverse na- 
ture. Observations indicate that 60 per cent are in spectro- 
scopic binaries (Milone & Latham 1992; Latham & Milone 
1996). Of those in binaries the orbits are a mixture of short- 
period and eccentric, long-period and eccentric, and long- 
period and circular. Blue stragglers are identified in a cluster 
colour-magnitude diagram (CMD) by their position above 
and blueward of the main-sequence turn-off. A popular ex- 
planation for their existence is that they are rejuvenated 
main-sequence (MS) stars that have gained hydrogen-rich 
material. Mass transfer in a binary from a companion to a 
MS star is an obvious means by which this can occur. There 
are three main scenarios (Kippenhahn, Wiegert & Hoffmeis- 
ter 1967) distinguished by the stellar type of the mass donor, 
MS star (Case A), sub-giant or red giant (Case B), or asymp- 
totic giant branch (AGB) star (Case C). Case A mass trans- 
fer can produce a BS in a very short-period Algol system, but 
in many cases ends in coalescence of the two MS stars when 
angular momentum is lost from the system and the orbit 
shrinks. This produces single BSs. Case B mass-transfer oc- 
curs in slightly wider binaries and results in a short-period 
circular binary containing a BS with a white dwarf com- 
panion. BSs in longer period binaries may be produced by 
Case C mass transfer or by accretion of material from an 
AGB star wind. The binary orbit will be circular except in 
some cases of wind accretion. However current models of 
binary tides and wind accretion only allow for accretion to 
occur and the orbit to retain an eccentricity in binaries with 
periods in excess of a few thousand days and the degree of 
accretion drops for wider orbits. This is a known problem in 
attempts to explain the existence of Barium stars in eccen- 
tric binaries (Pols et al. 2003). Thus binary evolution alone 
fails to explain the full range of BSs found in M67 and, in 
particular, the short-period eccentric binaries. It also fails 
to explain the so-called super-BS observed in M67 (Leonard 
1996) which has a mass greater than twice that of the cluster 
MS turn-off mass, Mto- Blue stragglers could also form as 
a result of direct collisions between MS stars although this 
is unlikely as the timescale for such an event to occur in an 
open cluster is too long (Press & Teukolsky 1977; Mardling 
& Aarseth 2001). 



1 M67 has been identified by the MODEST collaboration (Sills 
et al. 2003) as an ideal cluster for various groups work- 
ing on cluster evolution models to compare simulation tech- 
niques and use the observed data to calibrate the models: see 
|http : //manybody ■ org/modest/projects ■ html | 



In Hurley et al. (2001) we presented a semi-direct N- 
body model of M67. We showed that the combination of the 
cluster environment and a large population of binaries is 
able to generate formation paths for all of the BSs observed 
in M67. BSs were found in wide eccentric binaries because 
the BSs were exchanged into such orbits in three- and four- 
body encounters subsequent to their formation. Cases where 
two MS stars forming the inner binary of a triple system 
merged and remained bound to the third star were also ob- 
served and these produced BSs in short-period eccentric or- 
bits. Roughly half of the BSs formed during the simulation 
were the result of standard binary evolution in primordial 
binaries. The remainder owed their existence in some way to 
dynamical encounters. Some were the result of perturbations 
to primordial binaries. The orbits of these binaries are dy- 
namically altered so that mass-transfer begins whereas the 
stars would have remained detached if evolved in isolation. 
Collisions between MS stars at periastron in highly eccentric 
binaries also proved to be an efficient BS formation pathway. 
These eccentric binaries were the result of exchange interac- 
tions or perturbations to the orbits of primordial binaries. 
Super-BSs were also formed after the merger of three or more 
MS stars - an example of this sees a single BS formed via 
Case A mass transfer and then exchanged into an eccentric 
binary where it collides with its new MS companion. So this 
simulation was successful in showing that the cluster envi- 
ronment is able to boost BS production and that the diverse 
nature of the M67 BSs can be explained if they formed by a 
variety of mechanisms. However, the model did not generate 
enough BSs in binaries at any single moment - at most 25% 
of the BSs were in binaries when the model was at or near 
the age of M67. 

The main faihng of the Hurley et al. (2001) M67 model 
was that it is semi-direct. By this we mean that, owing 
to computational constraints, the model was not evolved 
by the direct A'^-body method for its entirety. At the time 
special-purpose hardware for computing gravitational forces 
was available in the form of the GRAPE-4 (Makino & Taiji 
1998) and while this had given a substantial increase in the 
performance of A^-body codes it still did not allow a com- 
plete model of M67 to be evolved from birth in a timely 
manner. The complication that M67 has a high frequency 
of binaries further reduced the simulation rate. The method 
employed by Hurley et al. (2001) was to take a population 
of 5 000 single stars and 5 000 binaries and evolve these from 
birth to an age of 2.5 Gyr according to a rapid stellar and 
binary evolution algorithm (Hurley, Tout & Pols 2002) - 
standard population synthesis. This evolved population was 
then used as input to the A^'-body code and evolved further 
to the age of M67 (about 4 Gyr, see next section) . The results 
of previous A-body simulations were used as a guide when 
selecting the stars and binaries for the model and to generate 
the density profile for the starting A'-body model. So every 
effort was made to account for the effect of processes such as 
mass-segregation during the 2.5 Gyr of cluster evolution that 
was skipped. However the situation was far from ideal. Since 
that time A'-body capabilities have improved markedly with 
the arrival of the next generation of special-purpose hard- 
ware, the GRAPE-6 (Makino et al. 2003), and its terafiops 
speed. In the meantime there have also been software and 
host-CPU speed-ups to aid with binaries (see Aarseth 2003). 
As a result binary-rich open clusters are now within reach 
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of direct A'^-body codes (Hurley & Shara 2003) as are larger 
single-star models (Baumgardt & Makino 2003). So we can 
now perform a direct model of M67. 

In this paper we attempt to create a direct A^-body 
model of the old open cluster M67. We discuss the model in 
comparison with observations of M67, looking at the overall 
structural properties and the nature of the stellar popula- 
tions. In Section|5|we describe the parameters of our starting 
model and in Sectional we look at the nature of the primor- 
dial binary population in some detail. We describe our sim- 
ulation method in Section |1| Section |K| contains the results 
of our first attempt at a model of M67 and then in Sectional 
we present our favoured M67 model. This is followed by a 
discussion and summary of our results in Section |7| 



2 THE INITIAL MODEL 

The goal is to evolve an A^-body model cluster from zero-age 
to obtain a model that resembles M67 at its current age in 
terms of both structural parameters and stellar populations. 
Unfortunately we cannot observe M67 at zero-age so in set- 
ting up our initial model we must be guided by observations 
of M67 at present in combination with the behaviour of pre- 
vious cluster models and also observations of young open 
clusters. In this section we discuss in turn our choices for 
the metallicity and age of M67, the binary fraction, the stel- 
lar initial mass function, the initial mass of the cluster, the 
tidal field and the density profile. The details of the primor- 
dial binary population are left for Section|21 We consider two 
models (blandly called Model 1 and Model 2) which differ 
only in the initial mass and the tidal field, and in the period 
distribution of the primordial binaries. All other choices are 
the same for both models. 

It is generally accepted that M67 stars are of solar com- 
position (Hobbs & Thorburn 1991; Tautvaisiene et al. 2000). 
There is no reason to believe that the metallicity of M67 
has varied over its life so we take Z — 0.02 for our initial 
model. The age of M67 is less settled with values quoted in 
the literature ranging from 3.2 ± 0.4 Gyr (Bonatto & Bica 
2003) to as high as 6.0 Gyr (Janes & Phelps 1994). In Hur- 
ley et al. (2001) we assumed an age of 4.2 Gyr based on 
fitting the isochrones of Pols et al. (1998) to the available 
data. Recently VandenBerg & Stetson (2004) derived an age 
of 4.0 Gyr for M67. We shall investigate behaviour around 
4 Gyr in the current study while bearing in mind that the 
cluster may be slightly younger or older. 

M67 is definitely a binary-rich cluster. Montgomery, 
Marschall & Janes (1993) found that at least 38 per cent of 
the cluster stars are binaries based on analysis of the pho- 
tometric binary sequence in the colour-magnitude diagram 
(CMD). A high binary frequency was confirmed by Fan et 
al. (1996) and Richer et al. (1998) who both found 50 per 
cent to be consistent with their data. How this translates to 
a primordial binary fraction depends somewhat on the pa- 
rameters of the primordial binaries and the evolution of the 
cluster. A good indicator is to look at the evolution of bi- 
nary fraction in previous open cluster models. Shara & Hur- 
ley (2002) performed simulations starting with 18 000 single 
stars and 2 000 binaries - a primordial binary fraction of 0.1. 
The orbital separations of the binaries were drawn from the 
distribution suggested by Eggleton, Fitchett & Tout (1989) 



with a peak at 30 au and limits of 0.03 au (about 6-R0) and 
30 000 au. This distribution was based on the Bright Star 
Catalogue (Hoffleit 1983) and is in agreement with the pe- 
riod distribution found by Duquennoy & Mayor (1991). Af- 
ter 4.5 Gyr of cluster evolution, Shara & Hurley (2002) found 
that the binary fraction was still close to 0.1. The evolution 
of binary fraction is shown in Figure for one of their mod- 
els. Also shown is the binary fraction for a simulation pre- 
sented by Hurley & Shara (2003) which began with 12 000 
single stars and 8 000 binaries. They used the same distribu- 
tion for the initial orbital separations but set the maximum 
at 200 au. This means the 40 per cent binaries would ac- 
tually represent a population of 55 per cent binaries if the 
upper limit were set at 30 000 au. In both cases the binary 
fraction remains similar as the cluster evolves. This can be 
understood because, while encounters may destroy soft bina- 
ries and a combination of dynamical hardening and binary 
evolution may destroy very hard binaries, the preferential 
escape of low-mass single stars from the cluster via tidal 
stripping counteracts these effects on the binary fraction. 
With a high fraction of primordial binaries residing in the 
high density core we do expect many binaries to be lost in 
binary-binary encounters. However, even in this region of the 
cluster the binary fraction does not drop because single stars 
are ejected from the core in the same encounters (Aarseth 
1996). The conclusion is that, for reasonable choices of the 
binary parameters, experience shows that in open cluster 
simulations the primordial binary fraction is well preserved 
as the cluster evolves. Therefore we assume a primordial bi- 
nary fraction of 0.5 and elaborate on the particulars of these 
binaries in the next section. 

Single star initial masses are chosen from the initial 
mass function (IMF) of Kroupa, Tout & Gilmore (1993, 
KTG93) between the mass limits of 0.1 and 5OA/0. All stars 
are assumed to be on the zero-age main sequence (ZAMS) 
when the simulation begins. This means that we neglect any 
age spread in the cluster stars that may be caused by dif- 
ferences in the pre-main sequence evolution timescales or by 
more than one epoch of star formation. We also assume that 
any residual gas from the star formation process has been 
removed. 

The number of stars in the initial model, or alterna- 
tively the initial mass of the cluster Mo, is an important 
quantity. To estimate it is difficult for us because determi- 
nation of the current mass of M67 from observations suff'ers 
from incompleteness and converting this to an initial mass 
requires knowledge of the mass-loss history from the cluster. 
A lower limit of approximately 8OOM0 for the current mass 
of M67 was set by Montgomery et al. (1993). Later Fan et al. 
(1996) estabfished that M67 has about lOOOM© in nuclear- 
burning stars with masses greater than O.5M0. We shall 
refer to this as the luminous mass, A^l. Hurley et al. (2001) 
showed that this corresponds to a total cluster mass, Al, of 
2 5OOAf0 and Mo ~ 3 5OOM0 if steUar and binary evolution 
are taken into account but the influence of the dynamical 
cluster environment is ignored. However, stars in a cluster 
are subject to mass segregation and this process, combined 
with the stripping of stars by the tidal field of the Galaxy, 
alters the mass function of the cluster over time to give a 
deficiency of low-mass stars. Looking at previous A'^-body 
data available to us (Shara & Hurley 2002; Hurley & Shara 
2003) we find that Ml = 1 OOOA^q is more likely to represent 
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a total mass of approximately 1 4OOM0 after 4 Gyr of stel- 
lar, binary and cluster evolution. We note that the limiting 
radius of the Fan et al. (1996) observations was 10 pc while 
the tidal radius of M67 is greater and therefore the actual 
mass of the cluster may be slightly higher. 

Next we must convert the current mass to an initial 
mass. Hurley et al. (2001) estimated that a starting model 
of 40 000 stars would be required to make a full model of 
M67. This was based on models with an escape rate pa- 
rameter fcc = 0.3 (see their equation 18). They also took 
the initial half-mass radius to be 1 pc which gave a starting 
model well within its tidal radius. A''-body models starting 
with 30 000 stars and no primordial binaries have recently 
been presented by Hurley et al. (2004). Analysis of these 
models shows that fee = 0.3 is once again suitable to de- 
scribe the average escape rate over the lifetime of the cluster. 
Using this value, and taking Mq — 15 OOOM0 for a cluster 
that fills its tidal radius initially and has a half-mass radius, 
Th, of 5pc (giving a half-mass relaxation timescale of about 
300 Myr), we find that approximately 2 OOOMq in stars is left 
after 4 Gyr. This neglects mass loss owing to stellar evolu- 
tion and there is also a degree of redundancy in the analysis 
but it shows that something of this order is required for the 
starting model. The 30 000 single-star models of Hurley et 
al. (2004) actually had M ^ 5 SOOM© at 4 Gyr which is too 
high for M67. However, it has been shown (Hurley & Shara 
2003) that the presence of a large number of primordial bina- 
ries increases the escape rate of stars and this is confirmed in 
our future paper on the general trends of binary-rich open 
cluster evolution (Hurley et al., in preparation). We note 
that the escape rate analysis of Hurley et al. (2001) does 
not consider the effect of binaries. Based on all this we start 
in our first attempt (Model 1) with 9 000 single stars and 
9 000 binaries which gives Mo — 14 4OOM0, comparable to 
that of the 30 000 single-star models. 

In Model 1 we assume for simplicity that the cluster 
is subject to a standard Galactic tidal field for the Solar 
neighborhood, a circular orbit with a speed of 220kms~^ at 
a distance of 8.5 kpc from the Galactic centre. For a mass of 
I44OOM0 this gives a tidal radius, rt, of 34 pc. M67 actually 
has a slightly eccentric orbit with a perigalacticon of 6.8 kpc 
and an apogalacticon of 9.1 kpc (Carraro & Chiosi 1994). We 
will therefore consider a second model (Model 2) starting 
with 12 000 single stars and 12 000 binaries on an orbit at 
8.0 kpc from the Galactic centre but leave further details of 
this to Section|Hl after we have presented the results for the 
first model. 

The density profile chosen for our starting model is that 
of a Plummer model (Plummer 1911; Aarseth, Henon & 
Wielen 1974). There is no reason to believe that a King 
model (King 1966), the major alternative, should describe 
the initial state of an open cluster as this set of models 
was originally conceived to describe dynamically evolved 
globular clusters. Kroupa, Aarseth & Hurley (2001) built 
a model of the Orion Nebula cluster which for simplicity 
used a Plummer density profile. They also assumed that the 
model contained twice the mass in gas as it had in stars and 
then after 0.6 Myr of A^'-body evolution the gas was removed 
smoothly. The result was not only a bound cluster but, after 
100 Myr, the projected radial density profile was a good fit to 
that of the Pleiades. In fact the choice is not likely to be cru- 
cial because previous A-body studies (e.g. Hurley & Shara 



2003) have shown that an initial Plummer model quickly 
evolves to resemble a King profile. The density profile in 
combination with the assumption that the stars are in virial 
equilibrium generates the initial positions and isotropic ve- 
locities of the stars. Formally the Plummer model extends 
to infinite radius so in the rare case of large distance (more 
than lOrh) a rejection is applied. The length-scale parameter 
is chosen such that the cluster fills its tidal radius from the 
beginning of the simulation. The position of the outermost 
star lies at a distance of 34 pc from the cluster centre. The 
choices for both starting models are summarized in Tabled 



3 PARAMETERS OF THE PRIMORDIAL 
BINARY POPULATION 

When setting up a population of primordial binaries we must 
choose how to distribute their orbital parameters. For each 
binary we need an orbital separation, a (or period, P), an 
eccentricity, e, and masses for the two stars (or one mass, 
which may be the total binary mass, Mb, and a mass ratio, 
q). The distribution function chosen for the orbital separa- 
tion or period is generally guided by observations with pop- 
ular choices being a flat distribution of log a (Abt 1983) or 
a log-normal distribution (Eggleton, Fitchett & Tout 1989; 
Duquennoy & Mayor 1991). For the eccentricity a choice is 
usually made between a thermal distribution (Heggie 1975), 
/(e) — 2e, or simply to have all orbits initially circular (al- 
though e = 0.0 is generally avoided for numerical reasons). 
When determining the masses of the component stars these 
may be uncorrelated, each mass is drawn independently from 
the same IMF, or a distribution of mass ratios may be as- 
sumed. 

In view of the relatively large number of blue stragglers 
found in M67, Hurley et al. (2001) performed a binary pop- 
ulation synthesis study to determine which combination of 
distribution functions for the orbital parameters of the pri- 
mordial binaries would maximize the number of blue strag- 
glers produced. Binaries were evolved according to the bi- 
nary star evolution (BSE, Hurley, Tout & Pols 2002) code. 
They found that choosing the orbital separation from a fiat 
distribution of log a in combination with correlated masses 
could possibly reproduce the M67 BS number provided that 
all BS producing binaries are retained by the cluster over 
its lifetime. So in this case standard evolution of primordial 
binaries is enough without the dynamical interactions be- 
tween cluster stars. However, the nature of the synthetic BS 
population does not correspond with that of the observed 
population. Only 25 per cent of the BSs are in binaries and 
all of these have circular orbits with almost all having an or- 
bital period of less than a year. Therefore, even if the cluster 
environment does not have a role to play in creating BS for- 
mation channels, it must be affecting the evolution of the 
binaries and the nature of the BS population. 

The case of uncorrelated component masses was ruled 
out by Hurley et al. (2001) because it leads to too few BSs 
by far. Their favoured population synthesis model had the 
Eggleton, Fitchett & Tout (1989) separation distribution 
with a peak at 10 an and a maximum of 200 au in com- 
bination with correlated masses and a thermal eccentricity 
distribution. The 200 au limit was based on the realisation 
that binaries wider than this would be weakly bound and 
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broken-up by dynamical encounters. This model predicted 
that about ten of the M67 BSs could be produced from stan- 
dard binary evolution of primordial binaries with 30% of the 
BSs in binaries, all with circular orbits. Hurley et al. (2001) 
then evolved their semi-direct model with the primordial bi- 
naries set up in this way. When the model reached an age of 
4.2 Gyr there were 22 BSs present. Unfortunately only one 
of the BSs present at 4.2 Gyr was still in a binary. 

A problem with the Hurley et al. (2001) semi-direct M67 
model was that no BSs were found in wide circular binaries 
at any time during the evolution while two are currently 
observed in M67 (Latham & Milone 1996). Exchange inter- 
actions do not tend to produce binaries in circular orbits so 
it would seem that the evolution of primordial binaries must 
produce any BSs in wide circular binaries. This can occur 
as a result of stable Case C mass transfer from an AGB star 
or via wind accretion. On inspection of the binary popula- 
tion in their starting model in combination with the setup of 
the BSE algorithm used to evolve the binaries, Hurley et al. 
(2001) discovered that no such binaries were expected. They 
found that if they repeated the binary population synthesis 
with the wind velocity reduced by a factor of two^ and the 
weaker criterion of Webbink (1988) used to detect the onset 
of dynamical mass transfer (as opposed to the default for the 
BSE algorithm) that the number of BSs predicted at 4.2 Gyr 
rose by 25 per cent for their favoured choice of distribution 
functions. More important the mix of the BS population 
was 22 per cent in wide circular binaries (up from per 
cent), 22 per cent in close circular binaries (down from 28 
per cent) and 56 per cent single (down from 72 per cent). 
We use a period of 1 000 d to make the distinction between 
close and wide binaries. In this paper we have used f3 — 1/8 
and the Webbink (1988) test for dynamical mass transfer 
in the BSE algorithm at all times, both when conducting 
population synthesis and within the A''-body code. 

An alternative approach to setting up the orbital pa- 
rameters of the primordial binaries is to use the method of 
pre-MS tidal evolution developed by Kroupa (1995b). Here 
it is assumed that the majority of low-mass stars form in 
aggregates that dissolve within 10 Myr (Kroupa 1995a) and 
that correlations between orbital parameters observed for 
short-period systems are the result of evolution in young 
binaries (less than 10^ yr) as the two accreting protostars 
interact. A birth period distribution can be constructed fol- 
lowing Kroupa (1995b) and tidal evolution is imposed on 
the eccentricity to replicate the observed deficiency of ec- 
centric orbits for short-period G dwarfs. A correlation be- 
tween the masses of the binary stars is also noted (Kroupa 
1995b). To investigate how this affects the production of BSs 
we set up the primordial binaries in the following manner. 
First the binary mass is chosen from the IMF of Kroupa, 
Tout & Gilmore (1991) and the component masses are set 
by choosing a mass ratio from a uniform distribution. This 
is the same method as used by Hurley et al. (2001) when 
they assumed correlated masses for the binary stars. Next 
a period is chosen according to the distribution given by 
Kroupa (1995b, equation 8) but with the generating func- 
tion given in Kroupa (1995a, equation lib) with S = 2.5, 



^ The wind velocity parameter, /3 was reduced from 1/2 to 1/8 
with the lower value more in keeping with observations 



r; = 45 and logPmin/d — 0. This gives logPmax/d — 8.43 for 
the distribution. An eccentricity is then chosen from a ther- 
mal distribution (Heggie 1975). If the periastron distance, 
Rpcii, for these parameters is less than five times the ZAMS 
radius of the primary, because pre-MS stars can typically 
be of this size (Kroupa 1995b), it is assumed that a collision 
would have occurred during the initial evolution and a new 
set of parameters is chosen. When this test is passed the 
birth eccentricity, Cb, is modified according to tidal evolu- 
tion by 

^ne, = ~(^y + ^ne^ {1) 

(Kroupa 1995b: equations 3b and 4), where A = 28, x = 0-75 
and e\ is the resulting eccentricity of the primordial binary. 

Generating up 500 000 binary stars and evolving them 
with the BSE algorithm we find that only 1.5 BSs per 10 000 
primordial binaries are predicted at 4 Gyr. This is in contrast 
to the model of Hurley et al. (2001) which predicts 11 BSs 
per 10 000 primordial binaries. Furthermore all of the tidally 
modified BSs are in wide circular binaries - no Case A or 
Case B mass transfer occurred. Recently Davies, Piotto & 
De Angeli (2004) have claimed that primordial binaries make 
mostly BSs in wide binaries and that the Case A/merger 
scenario is not dominant. This is based on the Preston & 
Sneden (2000) observations of field blue metal-poor (BMP) 
stars. However, from the Preston & Sneden (2000) data, it 
seems that if we assume that all BMP stars are BSs, then 
a possible mix of these BSs is 40 per cent from Case A 
evolution (single BS), 10 per cent from Case B (BS in short- 
period binary) and 50 per cent from Case C evolution (BS 
in wide binary). Furthermore, about 50 per cent of the wide 
binaries containing a BMP star have eccentric orbits and a 
possible explanation of these is that the BMP star (or BS) 
formed when the two stars comprising the inner binary of 
a triple merged. This would indicate that a large fraction 
of BMP stars were formed via the Case A/merger path. All 
we say on this for now is that much more could be done 
to constrain parameter distributions using population syn- 
thesis in combination with data from field stars and open 
cluster populations. We plan to look at this further in the 
near future. 

We will proceed with the Kroupa (1995b) tidal evolu- 
tion method, as described above, to generate the primordial 
binary population of our Model 1. This predicts about one 
BS from the primordial binaries so in effect the cluster envi- 
ronment must do all the work to create the observed popu- 
lation. Even though we expect that this will not be the case 
it is of interest to see if the dynamics can do it alone. 

In the previous section we mentioned that we will also 
evolve a second model which will start with more stars than 
Model 1 and orbit at a closer distance to the Galactic cen- 
tre (see Section |S| for details). For this model we also want 
the primordial binaries to be expected to produce BS num- 
bers of the order of what is observed in M67. The way we 
select the masses stays the same, Mb from Kroupa, Tout 
& Gilmore (1991) and n(g) = 1. However, instead of using 
the Kroupa (1995b) period distribution we simply select the 
initial separation from a flat distribution of log a and im- 
pose an upper cutoff at 50 au. For a binary with Mb = IM© 
this corresponds to a period of 1.3 x 10"" d. A population of 
50 per cent binaries chosen in this way actually represents 
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a population of 100 per cent binaries if all periods out to 
logPmax/d = 8.43 are allowed, as in Model 1. The assump- 
tion then is that all binaries born with a > 50 au quickly 
disassociate by weak encounters with other stars. Ignoring 
these non-interacting binaries of low binding energy will not 
affect the evolution of the cluster or the formation of stars 
such as BSs. At the lower end of the separation distribu- 
tion we reject any that is less than twice the sum of the 
ZAMS radii of the component stars. We then choose an ec- 
centricity from a thermal distribution (Heggie 1975), as be- 
fore, and modify it according to the Kroupa (1995b) tidal 
evolution test (Eq. if necessary. We note that Kroupa 
(1995a) claims that a fiat distribution of logP is consistent 
with pre-MS orbital data and that the Duquennoy & Mayor 
(1991) survey of binaries in the solar neighbourhood with 
G-dwarf primaries does not rule out a flat distribution of 
log a. Evolving 500 000 binaries set up in this way with the 
BSE algorithm gives 21 BSs at 4Gyr per 10 000 primordial 
binaries. The population synthesis predicts that 75 per cent 
of these should be single, 13 per cent in short-period circu- 
lar binaries, and the remaining 12 per cent in wide circular 
binaries. Thus the interest for Model 2 is in seeing if the 
combination of the cluster environment and the evolution of 
the primordial binaries can put the BSs in the required mix 
of living arrangements, while also maintaining the expected 
number of BSs. 



4 SIMULATION METHOD 

We use the NB0DY4 code (Aarseth 1999) to model the dy- 
namical evolution of star clusters. Basic integration of the 
equations of motion is performed by the Hermite scheme 
(Makino 1991) which employs a fourth-order force polyno- 
mial and exploits the fast evaluation of the force and its 
first time derivative by the GRAPE-6 (Makino et al. 2003). 
A time-step scheme comprising a series of hierarchical levels 
(McMillan 1986) allows each star to evolve on its own nat- 
ural dynamical timescale while forcing a block of particles 
to be advanced at each cycle so that efficiency does not suf- 
fer. Also, only one prediction is required for each block step. 
Regularization techniques are used to treat perturbed two- 
body motion in an accurate and efflcient manner (Mikkola 
& Aarseth 1998) with an extension to chain regularization 
(Mikkola & Aarseth 1993) to deal with compact subsystems 
of up to six bodies. A semi-analytical criterion developed by 
Mardling & Aarseth (2001) is utilized to detect and evolve 
stable hierarchical triple and quadruple systems which oth- 
erwise would prove extremely time consuming by direct inte- 
gration. Exchange interactions in encounters between single 
stars and binaries, or binary-binary encounters, where the 
member of a binary is displaced by an incoming star are 
included in this treatment. Direct collisions between stars 
(Kochanek 1992) and the formation of binaries in three- and 
four-body encounters are also allowed. Binary orbits may 
also become chaotic owing to perturbations from nearby 
stars and this scenario is modelled using the Mardling & 
Aarseth (2001) algorithm. 

In NBDDY4 stellar and binary evolution are performed in 
step with the dynamical evolution so that interaction be- 
tween these processes is modelled consistently (Hurley et 
al. 2001). Stellar evolution is included in the form of the 



single star evolution (SSE, Hurley, Pols & Tout 2000) algo- 
rithm. This is a package of analytical formulae fitted to the 
detailed models of Pols et al. (1998) that covers all phases 
of evolution from the ZAMS up to, and including, remnant 
phases. It is valid for masses in the range 0.1 — lOOM© and 
metallicity can be varied. The SSE package contains a pre- 
scription for mass loss by stellar winds that is utilized in 
NB0DY4. It follows the evolution of rotational angular mo- 
mentum for each star after the ZAMS spin orbital period 
is assigned according to a fit to the rotational velocity data 
of Lang (1992). All aspects of standard binary evolution are 
treated according to the BSE algorithm (Hurley, Tout & 
Pols 2002). Circularization of eccentric orbits and synchro- 
nization of stellar rotation with the orbital motion owing 
to tidal interaction is modelled in detail. Angular momen- 
tum loss mechanisms, such as gravitational radiation and 
magnetic braking, are also modelled. Wind accretion, where 
the secondary may accrete some of the material lost from 
the primary in a wind, is allowed with the necessary adjust- 
ments made to the orbital parameters in the event of any 
mass variations. Mass transfer also occurs if either star flUs 
its Roche lobe and may proceed on a nuclear, thermal or 
dynamical timescale. In the latter regime the radius of the 
primary increases in response to mass loss at a faster rate 
than the Roche lobe of the star. Stars with deep surface 
convection zones and degenerate stars are unstable to such 
dynamical timescale mass loss unless the mass ratio of the 
system is less than some critical value. The outcome is a 
common-envelope (CE) event if the primary is a giant. This 
results in merging or formation of a close binary, or a di- 
rect merging if the primary is a WD or low-mass MS star. 
On the other hand, mass-transfer on a nuclear or thermal 
timescale is assumed to be a steady process. Prescriptions to 
determine the type and rate of mass transfer, the response 
of the secondary to accretion and the outcome of any merger 
events are in place in BSE and the details can be found in 
Hurley, Tout & Pols (2002). 

One aspect that we should elaborate on here is the 
modelling of blue stragglers within BSE (and by associa- 
tion NB0DY4). If a MS star accepts mass from another star 
then it is rejuvenated. How this is done depends on whether 
the MS star has a radiative core (0.35 M/Mq 1-25), 
a convective core (Af > I.25M0) or is fully convective 
(Af < O.35Af0). As a MS star gains mass it evolves up 
along the MS to higher luminosity and effective tempera- 
ture. If the core is radiative the fraction of hydrogen that 
has been burnt is only slightly affected so that the effective 
age of the star decreases. In practice the age of the star is al- 
tered so that the fraction of its MS lifetime that has elapsed 
is unchanged by the change of mass. For MS stars with con- 
vective cores, and fully convective stars, the core grows and 
mixes in unburnt fuel as the star gains mass, so that the star 
appears even younger. The rejuvenation process is approx- 
imated by conserving the amount of burnt hydrogen and 
assuming that the core mass grows directly proportional to 
the mass of the star. We then take the remaining fraction 
of MS lifetime to be directly proportional to the remaining 
fraction of unburnt hydrogen at the centre to set the new 
effective age of the star. Owing to the increase in mass the 
remaining MS lifetime of the star has been reduced but it 
has been rejuvenated relative to other stars of its new mass. 

The merger or collision of two MS stars produces a new 
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MS star and we assume that the steUar material is com- 
pletely mixed in the process, with no mass lost from the 
system. The age of the new star is calculated on the as- 
sumption that, as stars evolve across the MS core hydrogen 
burning proceeds uniformly and that the end of the MS is 
reached when 10 per cent of the total hydrogen has been 
burnt. These stars are also rejuvenated relative to other 
stars of the same mass. The process of rejuvenation and 
subsequent evolution of the MS star are most likely ade- 
quate in the case of steady mass transfer although there is 
a lack of detailed studies to confirm this (Sills et al. 2003). 
The assumption that no mass is lost in a collision is in fair 
accordance with smoothed-particle hydrodynamics simula- 
tions of low- velocity collisions (Sills et al. 2001), but only 
a limited amount of mixing is found in these simulations so 
our assumption of complete mixing is questionable, and even 
more so in the case of slow mergers. In addition, the collision 
product is not initially in thermal equilibrium and requires 
a thermal timescale to contract to its MS state. During this 
time it has a higher probability to undergo additional inter- 
action (Lombard! et al. 2003). Furthermore the collision or 
merged product is, in most cases, rapidly rotating and needs 
to lose angular momentum in order to contract (Lombardi, 
Rasio & Shapiro 1996). These effects are not considered in 
our simplified model of post-collision BSs but would have an 
effect on the appearance of these stars. 

The great advantage of having the identical stellar and 
binary evolution algorithms in an A'^-body code and a popu- 
lation synthesis code is that we can evolve the same popula- 
tions inside and outside the cluster environment to quantify 
how the dynamical evolution affects the stellar populations. 
A'^-body simulations presented in this work were conducted 
on the 32-chip GRAPE-6 boards (Makino 2002) located at 
the American Museum of Natural History. 



5 RESULTS FOR MODEL 1 

The cluster of 9 000 single stars and 9 000 binaries set up in 
the manner described in Sections |5| and |3 was evolved from 
zero-age to an age of 5.8 Gyr when 1 000 stars remained. This 
simulation took approximately one month to complete. Fig- 
ure |5| shows the mass profile of the simulated cluster after 
4 Gyr of evolution. The mass remaining in the cluster at 
this time is 3175Af0. It is contained within a tidal radius 
of 21 pc. We note that the mass profile does continue be- 
yond the tidal radius because stars are not removed from 
the simulation until their distance from the density centre 
exceeds twice that of the tidal radius (cf. Terlevich 1987; 
Giersz & Heggie 1994). For this particular model the mass 
exterior to the tidal radius is 84M0, or 2.7% of the total 
mass. Also shown in Figure His the mass profile for the 
luminous mass which amounts to 1 987M0 . The luminous 
mass, that in stars above O.SMq burning nuclear fuel, con- 
tained within 10 pc is 1 73OM0 and it has a half-mass ra- 
dius of 3.0 pc. The half-mass radius for all stars is 4.9 pc. 
Fan et al. (1996) determined a luminous mass for M67 of 
1 OI6M0, rising to 1 27OM0 when corrected for binaries, for 
stars within 10 pc of the cluster centre. Our model has too 
much mass remaining to be a good representation of M67. 
It is not until 5.2 Gyr that we find Ml = 1 OOOM0 and we 



consider this too old to be relevant to M67. General results 
for this model are summarized in Table |5| 

At 4 Gyr the model contains only one blue straggler. 
This star resides in a non-primordial binary with an orbital 
period of 83 d and an eccentricity of 0.5. The mass of the BS 
is I.5M0 - compared to Mto = 1.32M0 - and the compan- 
ion is a O.64M0 MS star. Shghtly earher (3 950 Myr) there 
were two BSs and slightly later (4 050 Myr) there are three, 
so an average of two BSs at this age. The most BSs ob- 
served in the cluster at any one time, from an age of 3 Gyr 
onwards, is 5 at 4.5 Gyr. So the combination of the Kroupa 
(1995b) setup for the primordial binaries and the dynamical 
evolution of the cluster does not explain the actual M67 BS 
population. This is not a quirk of the model, either in terms 
of the initial conditions or statistics. We have performed a 
variety of simulations with the primordial binaries chosen in 
the same manner but have altered various characteristics of 
the starting model such as the King model for the density 
profile and the tidal radius filling factor. These models were 
evolved as part of a broader project to understand the evo- 
lution of binary-rich open clusters and will be described in 
detail in a forthcoming paper. In none of these models did 
we see substantial BS production. It has been noted previ- 
ously (Giersz & Heggie 1994) that statistical fluctuations in 
TV-body models may be reduced by averaging over the re- 
sults of many simulations. These fluctuations are amplified 
for small A^ and should not be of major concern in models 
of the size presented here (compared to the A'' — 500 mod- 
els discussed by Giersz & Heggie 1994). However, we have 
repeated Model 1 with an identical set up except for a dif- 
ferent seed for the random number generator. The evolution 
is similar - the mass remaining in the cluster after 4 Gyr 
is 2 957A/0 within a tidal radius of 20 pc and the luminous 
mass within the central 10 pc is 1 692M0 . At no time sub- 
sequent to an age of 3 Gyr are more than five BSs found in 
the model cluster. Therefore the nature of the primordial bi- 
naries in this type of simulation does not seem appropriate 
for a model of M67. Rather than continuing further with 
our analysis of Model 1 we instead turn our attention to 
Model 2. 



6 RESULTS FOR MODEL 2 - THE M67 
MODEL 

A fiaw of Model 1 was that it contained too much mass at 
4 Gyr to be considered a good model of M67. One way to 
reduce the mass remaining in a cluster at a certain age is 
to increase the strength of the tidal field in which the clus- 
ter evolves (Vesperini 1997; Baumgardt 2001). This reduces 
the tidal radius of the cluster which leads to an increase in 
the escape rate of stars. Associated with the reduction in 
mass is a lowering of the relaxation timescale (typically the 
half-mass radius will also be smaller) so that the cluster is 
dynamically older for a given physical age. The parameters 
given for the orbit of M67 in the Galaxy (Carraro & Chiosi 
1994) indicate an eccentricity of 0.14 and a time-averaged 
semi-major axis of 8kpc. So, when evolving Model 2, we 
will place the cluster on an orbit at 8.0 kpc from the Galac- 
tic centre. An orbital speed of 220kms~^ is still applicable 
at this distance (Chernoff & Weinberg 1990) so this remains 
unchanged and we keep the cluster on a circular orbit. 
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With the stronger tidal field we now expect the mass 
remaining after 4 Gyr of evolution to be less than that found 
for Model 1 and closer to the mass of M67, as desired. How- 
ever, we anticipate that the shorter relaxation time will drive 
the evolution faster than required to give the desired 25% 
reduction in cluster mass. So we make an associated increase 
in the particle number of the starting model to 12 000 single 
stars and 12 000 binaries. This gives Mo = 18 TOOA'/© within 
a tidal radius of 32 pc. 

The only other change in Model 2 compared to Model 1 
is in the set-up of the primordial binary population. As dis- 
cussed in Section|21a flat distribution of log a will be used to 
generate the orbital separations of the primordial binaries. 
This has a cutoff at 50 au which is comparable to the hard- 
soft binary limit (Heggie 1975) of about 40 au for Model 2 
with an initial half-mass radius of 3.9 pc. In all other re- 
spects the set up for Model 2 is the same as for Model 1. 
The two starting models are compared in Table Q From the 
12 000 primordial binaries in Model 2 we should expect 25 
BSs from binary evolution alone (see Section 1^. 

The mass profile for Model 2 at 4 Gyr is shown in Fig- 
ure|3 It has a total mass of 2 OSTA^q and a luminous mass of 
1 488M0 within a tidal radius of 15 pc. The luminous mass 
within 10 pc of the cluster centre (AfL.io) is now 1342AfQ 
which compares well with the binary corrected value found 
by Fan et al. (1996). Shghtly later in the evolution (4.1 Gyr) 
the model has A^l.io = 1 181Af0 with a tidal radius of 
14.5 pc. The half-mass radius of MS stars observed within 
10 pc was determined by Fan et al. (1996) to be 2.5 pc and 
our model has 2.7 pc for the same set of stars. So we now 
have a good match to the observed mass of M67 as well as 
the half-mass and tidal radii. General results for this model 
are summarized in Tabled 

In Figure |1] we once again show the mass profiles 
at 4 Gyr but this time normalized to the total mass of 
each profile. The luminous mass is more centrally concen- 
trated. Naively this result would be expected owing to mass- 
segregation when we consider that low-mass MS stars are ex- 
cluded from the luminous mass. However, the effect is coun- 
teracted somewhat by the added exclusion of white dwarfs 
(WDs) from the luminous mass because these are generally 
born in the central regions of the cluster and themselves 
are centrally concentrated relative to the entire population. 
Also shown in Figure |1] is the normalized mass profile of the 
BSs of which there are 20 in the model at 4 Gyr. These are 
more likely to be found in the centre of the cluster and have 
a half-mass radius of 1.1 pc. Fan et al. (1996) determined a 
half-mass radius of 1.6 pc for the M67 BSs. We elaborate fur- 
ther on the BSs and other stellar populations in Section r6.2l 

The cluster at 4 Gyr of age is a dynamically relaxed sys- 
tem - 13 half-mass relaxation times have elapsed in reaching 
this point. Figure|K|shows the behaviour of the number den- 
sity of stars contained within the core and the half-mass 
radius of the cluster as it evolved from to 4 Gyr. The 
core density of the starting model was 150 stars pc""^. This 
rose to a maximum of 330 stars pc~'^ after 3.5 Gyr. Here the 
core radius is the density-weighted value commonly used in 
theoretical models (Casertano & Hut 1985; Aarseth 2003) 
which is typically smaller than the core radius determined 
from observational techniques (Wilkinson et al. 2003). For 
binary-rich clusters there is no clear core collapse, at least 
not to the extent that we witness in simulations without pri- 



mordial binaries, where a high core density is required for 
binary formation (e.g. Hurley et al. 2004). Heating of the 
core by binaries occurs from the beginning in simulations 
with a large primordial binary population and this helps 
to keep the evolution of the core radius relatively regular 
and to avoid extreme fluctuations in central density. Stel- 
lar evolution mass-loss also leads to core expansion, espe- 
cially at early times when massive stars evolve away from the 
MS. The core density of the model at 4 Gyr is 83 stars pc~'^. 
The density of stars contained within the half-mass radius 
started at 50 stars pc~'^ and evolved to 4 stars pc~^ at 4 Gyr. 
Hence this was not a particularly dense system. At an age of 
5 Gyr the model cluster contains only 200 stars with a total 
mass of 26OAf0 and has almost reached the point of com- 
plete dissolution. In this paper the focus is on the results 
of the simulation at the age of M67 and the relevance of 
the model in improving our understanding of the observed 
properties of M67. We do not dwell on the details of the 
simulation in reaching an age of 4 Gyr. We shall present the 
long-term evolution of a binary-rich star cluster in another 
paper (Hurley et al., in preparation). 

6.1 Cluster Structure 

The coordinate system used in the A'^-body simulation has 
the cluster centre-of-mass as the origin, the X-axis directed 
away from the Galactic centre, the Y-axis in the direction 
of rotation about the Galactic centre and the Z-axis normal 
to the plane of the disk. In order to show how the model 
of M67 might appear on the sky we have rotated the model 
according to the Galactic coordinates of M67 {I = 215.68°, 
b — 31.93°: Bonatto & Bica 2003) so that the transformed 
YZ-plane becomes the observed plane (which we have also 
denoted the yz-plane). This view of the model is presented 
in Figure |S| . 

Bonatto & Bica (2005) used 2MASS data to construct 
a surface density profile for M67. This is shown in Figure |7| 
The limiting radius for the observations is 11.7 ± 0.6 pc and 
fitting a King (1966) surface density profile to the data gives 
a core radius of 1.14±0.13 pc within a tidal radius of 16±3 pc 
(Bonatto & Bica 2005). Using the same radial bins as for the 
observed M67 surface density profile, we also show the pro- 
file for the model in Figure Q The surface density of stars is 
greater in the model, especially in the central regions. How- 
ever, the observed profile includes only stars with J < 14.5 
which corresponds to stars with masses of about O.SMq or 
greater. Restricting the model profile to stars in this range, 
we find that agreement between the model and the observa- 
tions is better. The match is good in the 1 — 6 pc range which 
is the region that contains the half-mass radius of the cluster. 
Exterior to 6 pc the observed surface density drops to levels 
similar to that of the background counts (about Istarpc"'^) 
so it is not significant that the model surface density is lower 
in this region. We note that background contamination is nil 
in the A'^-body model. At the centre the model surface den- 
sity is still too high. The model is overdense by a factor of 
3 in the central bin. Fitting a King (1966) surface density 
profile to the model data gives a core radius of 0.6 pc so the 
model is clearly more centrally concentrated than observa- 
tions indicate for M67. Saturation effects could be lowering 
the observed star counts in the core of M67 but this is not ex- 
pected to be significant. We note that there is no noticeable 
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change in the model profile if we use different projections to 
construct it. 

A common method used in the analysis of observed clus- 
ter data is to construct a surface brightness profile. This 
can also be done with the A'^-body cluster data for which 
the mass, luminosity, stellar radius and position of each 
star is known, as well as whether or not the star is in a 
binary. First we calculate magnitudes and colours for each 
star using bolometric corrections provided by Kurucz (1992) 
and, in the case of WDs, Bergeron, Wesemael & Beauchamp 
(1995). Then the magnitude and position information is 
passed through the pipeline described by Mackey & Gilmore 
(2003) in their analysis of LMC clusters to give a projected 
surface density profile. The resulting data points for our M67 
model for the V-band in the XZ-plane are shown in Figure|S] 
A fit of the three-parameter Elson, Fall & Freeman (1987) 
model to the cluster profile is also shown and this gives a core 
radius of 0.66 pc. In this case all stars have been considered 
and the data is rather noisy. If we repeat the process but 
remove all stars with V < 12 in order to reduce saturation 
effects by bright stars and also remove all stars with V > 17 
to mimic a faint detection limit then we get the much cleaner 
profile shown in Figure |n| The fit to these data points gives 
a core radius of 0.64 pc. 

6.2 Stellar populations 

In Figure 1101 we present the colour-magnitude diagram 
(CMD) of Model 2 at 4.0 Gyr. There are 870 single stars 
and 1 325 binaries in the cluster at this time. The MS turn- 
off mass is 1.32M0. The cluster contains 2 968 MS stars, 57 
giants and subgiants, 2 naked helium (nHe) stars, 491 WDs 
and 2 neutron stars. Defining BSs as MS stars with mass in 
excess of the MS turn-off mass (by 2% or more) we find 20 
with masses in the range 1.4 — 2.1 Mq. These form the group 
of stars blueward of the the MS in the CMD. Nine of the 
BSs are in binaries. The BS population is discussed in detail 
in Sect. 7.2.1. There are six RS Canum Venaticorum (RS 
CVn) stars in the cluster. These are binaries with periods of 
1 ^ P/d ^ 14 that contain a cool subgiant star and a MS 
companion (Hall 1976). They are believed to be sources of 
X-ray emission. These and other expected X-ray sources in 
the cluster are discussed in Sect. 7.2.2. 

In the cluster at 4 Gyr are 226 single WDs. Most of these 
lie on the distinctive WD cooling track seen in the lower left 
corner of the CMD. There are 60 double- WD binaries, the 
majority of which are found in the CMD by their position 
above the sequence of standard single WDs. Of the double 
WDs, 28 have periods less than 1 d. The stars appearing 
in the region between the WD sequence and the MS are 
WD-MS star binaries in which the WD is still young and 
fairly hot. As the WD cools the binary moves across the 
CMD towards the MS as the MS star begins to dominate 
the appearance of the binary. Further particulars of the WD 
population are discussed in Sect. 7.2.3. There is also one 
cataclysmic variable (CV) in the cluster at this time, located 
atV = 23.06 and (B-V) = 0.53 in the CMD. It comprises a 
O.IM0 MS donor star in a 0.55 d circular orbit with a O.3M0 
helium WD. The CV phase began at an age of 1.9 Gyr when 
the MS star mass was O.I8M0. 

A summary of selected stellar population results for 
the M67 model is given in Table |3| along with results from 



Model 1. There are two extremely blue stars at « 12 in the 
CMD. These are not BSs but binaries comprised of evolved 
nHe stars with MS star companions. Both evolved from pri- 
mordial binaries in which the nHe star was produced in a 
CE event initiated by the progenitor of the nHe star filling 
its Roche lobe while on the early AGB. However, the fainter 
of the two binaries would not exist in this state without 
having experienced a significant perturbation to its orbit. 
It originated as a primordial binary with an eccentricity of 
0.5 and a period of 30 902 d. Evolved in isolation the two 
stars would not have become close enough to interact and 
the result would be a WD-MS star binary with an orbital 
period of about 300 yr at 4 Gyr. Within the cluster environ- 
ment the binary received a perturbation to its orbit from 
a third star, in a fiyby encounter after 350 Myr, while still 
a MS-MS star binary. This resulted in a slight decrease in 
orbital period but more important pumped the eccentricity 
up to 0.95 and subsequent tidal evolution brought the stars 
close enough for mass-transfer to begin. The CE event then 
reduced the orbital period further so that the nHe-MS star 
binary observed at 4 Gyr has P — 12.6 d and is circular. The 
other nHe-MS star binary has a period of 3.5 d. 

Another anomalous star in the CMD lies towards the 
base of the giant branch (GB) but below the subgiant branch 
atV = 13.32 and {B - V) = 0.83. This is a single star with 
a mass of I.9IM0 recently created in a CE merger event. 
At an age of 3.8 Gyr a 1.37M0 subgiant primary filled its 
Roche lobe and began steady Case B mass transfer on to 
its O.83M0 MS star companion. When the primary reached 
the end of the subgiant phase (about 50 Myr later) its mass 
had been reduced to O.98M0 and the companion mass was 
I.22M0 owing to the mass-transfer being conservative. At 
this stage the envelope of the primary was fully convective 
and the entire envelope overfiowed the Roche lobe on a dy- 
namical timescale to create a common envelope. The mass 
ratio inversion had not been enough to avoid this. The gi- 
ant core and the MS star spiralled together and expelled 
O.29M0 of the envelope via dynamical friction before merg- 
ing to form a I.9IM0 giant with a core mass of O.I6M0. This 
star is more massive than the normal giants in the cluster 
at this time but has a core mass less than expected for a 
giant of this mass and in fact less than the core mass of the 
stars residing at the base of the cluster GB. The reduced 
core mass is the result of the mass loss experienced by the 
subgiant progenitor which restricted its core mass growth 
as it evolved across the subgiant branch. This is the cause 
of the merged giant lying below the subgiant branch. The 
increased mass of the giant will cause it to remain on the 
blue side of the standard cluster GB as it continues its giant 
evolution. 

6.2.1 Blue stragglers 

Table m provides a list of the 20 BSs in the M67 model. All 
of the eleven single BSs were produced from the merger of 
two MS stars in a primordial binary. In nine of these cases 
BS formation was via the onset of Case A mass-transfer and 
the eventual coalescence of the stars as angular momentum 
was removed from the system. For seven of these the evolu- 
tion proceeded as if the binaries were evolved in isolation. 
The cluster environment did not interfere. In one case (star 
#2411), a perturbation to the orbit hastened the onset of 
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Case A mass-transfer so that a BS that would not have been 
created until after 4 Gyr was formed earlier and observed at 
4 Gyr. In another case (#3021) mass transfer was delayed by 
the involvement of the binary in a temporary exchange in- 
teraction which increased the orbital period. This BS would 
not have been observed at 4 Gyr without the interference of 
the third star as it would have already evolved to become 
a WD. On the other hand, if the exchange interaction had 
caused a larger period increase the BS may have formed af- 
ter 4 Gyr or not at all. The remaining two of the single BSs 
formed from initially long-period binaries that had their ec- 
centricity pumped up to 0.99 by weak fly by encounters. This 
led to a collision of the MS star components at periastron 
and creation of a BS. In neither binary would the compo- 
nent stars have become close enough to interact and merge 
without the intervention of a third star. Descriptions of the 
evolution pathways for the four single BSs that were affected 
in some way by the cluster environment are given in Table |S] 

Descriptions for the nine binary BSs are also given in 
Table |^ Two of these are primordial binaries but both 
have been affected by interactions with other cluster mem- 
bers. This explains the eccentric nature of the orbits. In 
one case the BS was formed in a wide circular binary as 
a result of Case C mass-transfer and a subsequent pertur- 
bation induced the eccentricity observed in the orbit. The 
other primordial binary became involved in a four-body in- 
teraction and one of the binary components collided with 
another member of the subsystem to form the BS which 
then emerged bound to its original companion. The remain- 
ing BS-binaries are non-primordial and their formation in- 
volved an exchange interaction at some point. In four of 
these cases a primordial binary became part of a three- or 
four-body system and perturbations to the orbit drove the 
eccentricity of the binary towards unity so that the stars 
collided. The merged star, a BS, was then exchanged into a 
binary. One of these cases then underwent a second collision 
in an eccentric binary while another was subsequently per- 
turbed. A BS binary was formed from Case-C mass transfer 
in a primordial binary but the BS was then exchanged into 
a new binary and, just prior to 4 Gyr, Case B mass-transfer 
began in this binary. This further increased the mass of the 
BS. Another case saw an exchange interaction form a binary 
which evolved to a state of Case-A mass transfer followed 
by coalescence to a BS which was later exchanged into a 
new binary. The final case involved an exchange interaction 
followed by two collision events. So the binary BSs were 
formed by a variety of means and this resulted in a mix of 
orbital parameters: short-period and circular, short-period 
(less than 1 000 d) and eccentric, long-period (greater than 
lOOOd) and eccentric (see Table 2] for full details). 

Only seven of the 20 BSs in the model evolved from 
unperturbed primordial binaries. Of the remaining 13, eight 
BSs were formed by collisions that were induced by three- 
or four-body interactions, or by perturbations that drove up 
the eccentricity to almost unity. These eight could not have 
formed outside the cluster environment and the same is true 
for the case A merged star that formed after an exchange. 
The other four would have become BSs by binary interac- 
tion alone but two of these would not have been observed 
as BSs at an age of 4 Gyr. Hence, for approximately half 
the BS population the dynamical cluster environment was 
instrumental in producing them. 



As mentioned earlier the half-mass radius for the BSs 
in the M67 model is 1.1 pc. This is much less than the half- 
mass radius for all stars (3.8 pc) and indicates that the BSs 
are a centrally concentrated population which is in agree- 
ment with their observed distribution in M67 (Fan et al. 
1996). The fact that BSs are more likely to be found in the 
centre of the cluster than the outer regions is conflrmed by 
Figure [TTI which splits the CMD into two regions, stars con- 
tained within the half-mass radius and stars exterior to this. 
We see that only two of the BSs in the model at 4 Gyr are 
found outside of the half-mass radius and one of these was 
situated well within the half-mass radius when it formed"^. 
The reasons for this are twofold. Primordial binaries which 
are to produce a BS must have a mass in excess of the cur- 
rent MS turn-off mass (1.3 A/0) which itself is greater than 
the average mass of the cluster stars (O.6M0). So the pro- 
cess of mass-segregation acts to concentrate these binaries at 
the centre of the cluster and the BSs themselves, when cre- 
ated, also tend to sink towards the centre. We also have BSs 
created in binaries formed from exchange interactions and 
these interactions are more likely in regions of high density 
(Heggie, Hut & McMillan 1996) which favours the centre of 
the cluster (see Figure (KJ. 

The 20 BSs in the model at 4 Gyr is less than the 25 
predicted when evolving the primordial binaries in isolation 
(see Section|^. Considering that only 11 of the 20 BSs would 
have formed without help from the cluster environment, this 
means that 14 potential BSs at 4 Gyr were lost during the 
simulation. So in addition to creating BSs the cluster envi- 
ronment is just as productive in destroying potential BSs, 
maybe even more so. The majority of lost BSs were the result 
of perturbations to primordial binaries that hardened the or- 
bit and brought forward the onset of Case-A mass transfer. 
BS formation still occurred in these binaries but much earlier 
than it would have otherwise and this caused the BS phase 
to have ended prior to 4 Gyr. This is compensated to some 
extent by somewhat wider binaries being hardened into or- 
bits that allow them to experience Case-A mass transfer and 
coalesce but otherwise would not have done so within 4 Gyr, 
and by interactions that delay the onset of mass transfer as 
in the case of #3021. Other possibilities include an exchange 
interaction destroying the primordial binary or the binary 
being ejected from the cluster. Both of these events are rare 
compared to the hardening scenario. The former because the 
pre-exchange binary is short-period and the latter because 
the binary is relatively heavy. 

There are no BSs in long-period (nearly) circular bina- 
ries present at 4 Gyr whereas there are two observed in M67: 
P = 1221d, e = 0.09 and P = 1154d, e = 0.07 (Latham 
& Milone 1996). However, BSs in such binaries were present 
in the simulation at other times. An example is a 1.8Mq 
BS in a circular orbit of period 1 445 d with a WD compan- 
ion at 2.6 Gyr. The proto-BS initially accreted mass from 

^ This BS lies outside of the tidal radius at 4 Gyr, 19.4 pc from the 
cluster centre. It is counted as a cluster member because its orbit 
is such that it remains bound to the cluster and subsequently 
returns within the tidal radius. The particulars of this BS are 
provided in Table |3 (BS #2565) and its evolution pathway is 
detailed in Table |^ The radial position of this BS explains why 
the normalized mass profile of the BSs in Figure|^does not reach 
unity within the tidal radius. 



N -body model of M67 11 



the stellar wind of its AGB-star companion and then grew 
even more in mass with the onset of Case-C mass transfer, 
with dynamical timescale mass-transfer and CE evolution 
avoided because the stellar wind had significantly decreased 
the mass of the AGB star by this stage. Mass transfer ceased 
when the AGB star exhausted its envelope and became a 
WD. Also, BS #1378 observed at 4Gyr in a long-period ec- 
centric binary was originally in a wide circular binary when 
formed at T = 2 Gyr. It would have been found in this state 
at 4 Gyr if the orbit had not been perturbed by a passing 
star at 2.5 Gyr. So in a sense it is simply bad luck that we 
did not observe any BSs in long-period circular binaries after 
4 Gyr of evolution. The same can be said for the incidence 
of super-BSs. M67 is observed to have a super-BS of mass 
3Mq (Leonard 1996) but our model at 4 Gyr does not con- 
tain any BSs with mass greater than twice the MS turn-off 
mass of the cluster. Super-BSs do form in the simulated clus- 
ter, a total of five for the entire simulation including one at 
3.9 Gyr when the model contained 22 BSs. This super-BS of 
mass 3.2M0 was found in a binary of period 4.6 d and ec- 
centricity 0.45 with a I.2M0 MS star companion. The star 
began its life as a I.2M0 MS star with a I.3A/0 compan- 
ion in a circular primordial binary with P = 1.2 d. After 
2.7 Gyr the more massive MS star filled its Roche lobe and 
Case-A mass transfer proceeded until 3.1 Gyr when the stars 
coalesced to form a 2.5M0 single BS. The BS was situated 
well within the core of the cluster (0.1 pc from the centre) 
and at 3.4 Gyr exchanged itself into an existing binary. Its 
new companion was a O.7M0 MS star and the orbit had an 
eccentricity of 0.1 and a period of 1.3 d. The orbit quickly 
circularized owing to tidal forces and Case-A mass transfer 
began at 3.6 Gyr with coalescence almost immediately after- 
wards. This newly formed 3.2M0 super-BS was exchanged 
into a wide binary with a I.3M0 companion at 3.7 Gyr and 
this binary was then involved in a short-lived four-body sys- 
tem that created the short-period eccentric binary observed 
at 3.9 Gyr. The super-BS evolved off the MS about 50 Myr 
later and just missed being included in the M67 model. This 
evolution example also serves to highlight how BSs may be 
formed in short-period eccentric binaries. M67 is observed to 
contain a BS in an orbit with period 4.18 d and eccentricity 
of 0.2 (Milone & Latham 1992). 

The orbital parameters of all BS-binaries created in 
Model 2 at an age of 2 Gyr or later are shown in Figure 1121 
Also shown are the six M67 BS-binaries with known orbital 
solutions. Figure ITHl shows the total number of BSs and BSs 
within binaries as a function of cluster age up to 4.2 Gyr. We 
note that the numbers we have focussed on at 4 Gyr are typi- 
cal of the cluster over the preceeding Gyr or so. The increase 
in BS number seen after 2 Gyr corresponds to an increase 
in central density as the cluster becomes dynamically more 
evolved. The BSs found in binaries during this timeframe 
are more likely to be in non-primordial binaries which is the 
opposite of what is found earlier in the evolution. In Fig- 
ure^] we also show the number of BSs in circular binaries 
with P < 100 d. Only one BS binary of this type is found in 
the cluster at 4 Gyr but at earlier ages these binaries, formed 
primarily from Case B mass transfer, were dominant. Of the 
BS binaries formed prior to an age of 2 Gyr 82 per cent were 
circular with an orbital period of 100 d or less which con- 
trasts with 23 per cent of the BS binaries in Figure lT^ being 
of this type. The number of these binaries declines with age 



and this is linked to the destruction of primordial binaries 
as the cluster ages. 

We have seen that formation scenarios for all of the 
various BSs observed in M67 exist in Model 2. The model 
at 4 Gyr has a binary fraction of about 0.5 for the BSs which 
is close to the fraction of 0.6 observed. Also, the ratio of BSs 
to the number of MS stars within two magnitudes of the MS 
turn-off (13 < V < 15) is 0.18 which matches well with the 
observed value of 0.14 (Ahumada & Lapasset 1995). In fact 
this indicates that the model is over-producing BSs although 
if we look at raw numbers the model has 20 BSs at 4 Gyr 
compared to 28 for M67 (Hurley et al. 2001). However, if we 
look at the CMD of M67 (e.g. Fig. 2 of Hurley et al. 2001; 
Montgomery, Marschall & Janes 1993) the BSs appear to 
form two distinct groups, an obvious group of 10 BSs with 
V < 12 and the remainder much closer to the MS turn-off 
position. Some of these less obvious BSs actually lie below 
the MS turn-off by almost a magnitude. The observed BSs 
are identified by their position in the CMD whereas in the A*'- 
body model we have much more information about the stars 
and have used mass as the determining factor. Inspection of 
Figure nm and comparison with the positions of the observed 
M67 BSs reveals that there are at least five stars near the 
MS turn-off of the model that we may have classified as BSs 
if using CMD position as the determining factor. Counting 
these stars as BSs would give good agreement of raw BS 
numbers for the model and M67. Also, two of the M67 BSs 
have proper motion membership probabilities of less than 
80 per cent (Girard et al. 1989). 

We now focus on a comparison with the group of ten 
M67 BSs with V < 12. This sample has been well studied 
by Milone & Latham (1992) revealing that six are single 
and four are in binaries. The orbital parameters of the BS- 
binaries are P = 4.2d and e = 0.2, P = 846 d and e = 0.5, 
P = 1003 d and e = 0.3, and P = 1 221 d and e = 0.1. 
The simulated cluster at 4 Gyr also has a pronounced group 
of 15 BSs (those with \B - V\ < 0.45 in Figure ITIl. Eight 
of these are single and seven are in binaries. If we further 
restrict this sample to V < 12 then we have eleven BSs with 
six in binaries. The widest of these six (#1613 and #3835) 
may not be detected as having a companion if observed so 
the mix of BSs in this sample could easily switch to 7 single 
BSs and 4 binary BSs. Either way, considering the stochastic 
nature of BS formation, the number of BSs in the restricted 
sample is a good match to that of the observed group of BSs 
as is the ratio of single BSs to binary BSs and the orbital 
combinations of the binary BSs. 

6.2.2 X-ray sources 

X-ray observations of old open clusters (van den Berg et al. 
2004) and globular clusters (Pooley et al. 2003) have proved 
to be very efficient at detecting short-period interacting bi- 
naries. For open clusters there is a well-defined connection 
between age and X-ray activity (Randich 1997) where the 
latter declines with age owing to the spin-down of late-type 
stars which are initially rapid rotators. Therefore, X-ray de- 
tections in a cluster such as M67 are expected to be associ- 
ated with stars that have been spun up in some way, typi- 
cally by residing in a short-period binary where the spin pe- 
riod of the star is kept synchronized with the orbital period 
by tidal forces. An example of such a binary is an RS CVn 
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system which contains a subgiant and a MS star companion 
in a short-period orbit. BY Draconis binaries are chromo- 
sphericaUy active systems containing a late-type (spectral 
type F or later) MS star primary with a MS star companion 
and may also explain X-ray activity in old populations. 

A Chandra observation of M67 was recently reported by 
van den Berg et al. (2004). They detected 158 X-ray sources 
(with a limiting flux corresponding to an X-ray luminosity 
of Lx ~ 10^* ergs~^). Optical counterparts that are proper- 
motion members of M67 were found for 25 sources and a fur- 
ther 12 sources had optical counterparts that are believed to 
be cluster members based on their position near the lower 
end of the M67 main sequence. Ten of the proper-motion 
members are binaries with periods less than 12 d and con- 
tain subgiant or MS stars. Approximately six of these appear 
to be classical RS CVn binaries of which three were detected 
earlier by Belloni , Verbunt & Mathieu (1998). The twelve 
sources near the lower end of the MS are possible BY Dra- 
conis binaries. 

Table |S| lists the parameters of the six RS CVn bina- 
ries observed in our M67 model at 4Gyr. Five of these are 
the result of standard primordial binary evolution while one 
(#1568) was formed in an exchange interaction. The number 
is a good match to the M67 observations although possibly 
we should not count the system in which the primary is fill- 
ing its Roche lobe (#2633) or the system with a period of 
20 d (#2383) when making the comparison as these fall out- 
side of the classical RS CVn period range as defined by Hall 
(1976). 

To investigate the incidence of BY Draconis binaries 
we look at all MS-MS binaries where the primary is F-type 
or later (Mi ^ l.OM©) and the orbital period is 12 d or 
less. This cutoff in period is used in observational work (e.g. 
van den Berg et al. 2004) and is thought to ensure that all 
binaries in the sample are synchronously rotating at the age 
of M67. Using the relation between X-ray luminosity and 
rotation rate, Q,, proposed by Walter (1982), 

log Lx/iboi = -3.14-0.16/0, (2) 

where Lboi is the bolometric luminosity of a star and Q, is in 
units of rotations/day, we can estimate an X-ray luminosity 
for each of the BY Draconis binaries. These are plotted in 
Figure 1141 as a function of orbital period. The upper limit 
to the X-ray luminosity at any period is given by a I.OMq 
star in synchronous rotation with its orbit - to emphasize 
this we have plotted the behaviour of Equation 2 for such 
a situation (upper solid line in Figure [T^ . Correspondingly 
the lower limit occurs for a O.IMq star - the lowest mass of 
star that we considered (see the lower solid line in Figure lT^ . 
A small number of systems in our M67 model have Lx less 
than this lower limit and these are binaries where both stars 
are less massive than O.2M0 and synchronous rotation is 
yet to be achieved. We have also included data points for 
the BY Draconis systems identified by van den Berg et al. 
(2004) for which orbital periods are known (four of twelve). 
We find 172 MS-MS systems with Lx > 10^*ergs"\ Also 
shown in Figure 1141 are 19 short-period MS- WD binaries 
which have X-ray luminosities in this range and for which 
the WD component has cooled sufficiently that the binary 
appears near the MS in the CMD. 

We have an over-abundance of potential MS X-ray 
sources compared to the Chandra detections for M67. How- 



ever, the Chandra pointing involved six CCD detectors with 
each having an area of 8.4' x 8.4' so that the coverage was 
only about l/30th of the full area of M67 (for a tidal radius 
of 60'). Also, the sensitivity of the Chandra observations 
decreases towards the edge of the detector and the limit- 
ing luminosity will be greater than lO'^* ergs~^ for the outer 
detectors (M. van den Berg, private communication). The 
centre of the cluster was covered by Chandra so we do not 
expect the RS CVn sample to suffer greatly from incom- 
pleteness as these are relatively massive binaries but BY 
Draconis binaries will not necessarily reside in the central 
regions of the cluster, especially those with cool primaries. 
Therefore there may be a substantial population of these 
systems yet to be observed. 

6.2.3 White dwarfs 

Hurley & Shara (2003) used A^-body models to conduct a de- 
tailed investigation of the behaviour and appearance of the 
white dwarf population in dense star clusters. They found 
that the mass fraction of WDs can be significantly enhanced 
by the dynamical evolution of the cluster - by as much as a 
factor of 2 in old open clusters. This enhancement is relative 
to the mass fraction expected if the same populations were 
evolved with basic population synthesis and with the same 
stellar and binary evolution algorithms used in the A-body 
code. A combination of mass segregation and the presence 
of a tidal field mean that, as a cluster evolves, low-mass MS 
stars are preferentially stripped from the cluster and the 
WD mass fraction diverges from that of the field popula- 
tion. We find the same behaviour in our M67 model. The 
expected WD mass fraction at 4Gyr is /wd =0.1 when 
we perform population synthesis on the initial population. 
However, the M67 A-body model has /wd = 0.15 at 4Gyr. 
The average mass of the cluster WDs starts high and then 
decreases with time as progressively less-massive WDs are 
born. After 1 Gyr the WD average mass is O.82M0 and after 
4Gyr it is 0.64Mq. On the other hand the average mass of 
MS stars is O.52M0 for the starting model and initially de- 
creases owing to evolution of the massive MS stars followed 
by an increase to O.56M0 at 4 Gyr because of tidal strip- 
ping. In other words, the average mass of MS stars is nearly 
constant. Based on these values we would expect the WDs 
to be preferentially found near the centre of the cluster as a 
result of mass segregation. This expectation is strengthened 
when we consider that WDs are born from stars that were 
previously the most massive MS stars in the cluster. What 
we observe in our M67 model is that the half-mass radius of 
the WDs is 0.64 pc which is indeed centrally concentrated 
when compared to a half-mass radius of 3.8 pc for all stars. 
If we look at only single WDs the half-mass radius increases 
to 1.28 pc (for comparison the half- mass radius of the cluster 
binaries is 0.37 pc). The upshot of all this is that low-mass 
MS stars are more likely to reside in the outer regions of 
the cluster and are thus more vulnerable to escape from the 
cluster than WDs. 

Figure 1151 shows the variation of WD mass fraction as 
a function of radial position in the M67 model. We see that 
the value of /wd depends strongly on which portion of the 
cluster is observed. Somewhat surprisingly the greatest en- 
hancement compared to the population synthesis value is 
found in the 6 to 10 pc region, exterior to the cluster half- 
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mass radius. There is less but noticeable enhancement in the 
core which should contain a significant population of mas- 
sive old WDs as well as double WDs and less-massive young 
WDs. However, the central regions are also dominated by 
massive stars and binaries of all types which reduces /wd- 
This also means that young WDs move outwards from the 
centre as they cool. Outside the half-mass radius there is a 
lack of massive stars and binaries, as well as young WDs (see 
also Figure 111^ . This increases the relative mass of the WD 
population. There is one radial bin near the half-mass ra- 
dius where /wd is similar to the population synthesis value. 
Unfortunately the behaviour of /wd in the vicinity of this 
region is rather erratic and this makes it difficult to suggest 
that observations to determine /wd should be conducted 
near the half-mass radius. In fact our results indicate that 
M67 is dynamically old enough that a measurement of /wd 
cannot be used to yield information about the IMF. Also 
shown in Figure 1151 is the value of /wd calculated if we 
ignore WDs in binaries with non-WD companions. So this 
value corresponds to single WDs and double WDs found on 
or near the WD sequence in the cluster CMD which are the 
WDs most likely to be observed in a real cluster. 

Richer et al. (1998) derive a WD mass fraction of 0.09 
for M67. This is based on finding 85 WDs down to the termi- 
nation point of the WD cooling sequence (their deep survey 
reached V = 25). Correcting this number for WDs hidden 
in binaries and taking a 50 per cent binary fraction they 
estimate that there could be as many as 150 WDs in M67. 
Assuming an average mass of 0.7Mq for the WDs and using 
a calculated total mass of 1 O8OAf0 for M67 they arrive at 
the quoted value of /wd. Richer et al. (1998) note that this 
number appears to be low when compared to the number 
of giants they observed. Based on a population of 87 giants 
they expected about 60 WDs with a cooling age less than 
1 Gyr but found only 24. Our M67 model contains 226 single 
WDs, 60 double- WD binaries and 145 WDs in binaries with 
a non-WD companion. Of the single WDs 54 have a cooling 
age less than 1 Gyr. There are 22 double WDs with at least 
one bright component and 33 bright WDs contained in other 
binaries (some of these with low-mass MS companions may 
appear near the WD sequence). The model also contains 59 
giants which gives an approximate 1:1 ratio of bright WDs 
to giants if we consider only the single WDs. If our model 
is to be believed this suggests that the observations of WDs 
in M67 are incomplete. 

Nine double- WD binaries with a total mass in excess 
of the Chandrasekhar mass (Mch = 1.44M0) and merger 
timescales owing to gravitational radiation of less than the 
age of the Galaxy (about 13 Gyr) are produced in Model 2. 
One of these contains two oxygen-neon (ONe) WDs and 
three have an ONe and a carbon-oxygen (CO) WD compo- 
nent. The outcome of WD- WD mergers involving an ONe 
WD is believed to be an accretion-induced collapse (AIC: 
Nomoto & Kondo 1991) to form a neutron star remnant. 
This is the assumed outcome when such an event occurs in 
an NB0DY4 simulation. Mergers of CO-CO WD binaries with 
M > A'/ch are possible causes of Type la supernovae (Yun- 
gelson & Livio 2000) but may instead lead to AIC and the 
formation of a neutron star (Saio & Nomoto 1998). Shara 
& Hurley (2002) found that the incidence of these possible 
Type la events increased by as much as a factor of 10 in 
the environment of an open cluster. This result was based 



on simulations of 20 000 stars with a 10 per cent primordial 
binary population and the Eggleton, Fitchett & Tout (1989) 
orbital separation distribution used to set up the binaries. 
From the 12 000 primordial binaries in our simulation, pop- 
ulation synthesis with the BSE code predicts that we should 
get three super-Chandrasekhar CO-CO WD mergers. So the 
five that we did find represents an enhancement but not of 
the order found by Shara & Hurley (2002). However, it is too 
early to tell if the degree of enhancement depends on pri- 
mordial binary fraction. The two additional CO-CO mergers 
in Model 2 are the result of perturbations to the orbits of 
primordial binaries that brought the component stars closer 
together. In each case a CO-CO WD binary would have 
formed regardless of the perturbation but the orbital pe- 
riod would have been greater and merging would not have 
occurred. 

At 4 Gyr one of the ONe-CO WD merger candidates 
remains in the cluster. It is not expected to merge for an- 
other 8 Gyr which is after the cluster will have completely 
dissolved. The other two ONe-CO WD binaries merged and 
formed neutron stars at ages of 63 and 75 Myr. The ONe- 
ONe WD binary formed after 61 Myr and the components 
were expected to merge shortly afterwards but in the interim 
the binary was ejected from the cluster. The five CO-CO 
WD binaries merged prior to 4 Gyr, the first after 110 Myr 
of evolution and the last at 2 666 Myr. Also formed during 
the simulation were ten double helium WD binaries with 
merging timescales less than 13 Gyr. Four of these escaped 
from the cluster prior to 4 Gyr. The other six remain in 
the cluster with periods of 0.04 to 0.08 d and three are ex- 
periencing steady mass transfer. In the BSE and NB0DY4 
codes mass transfer from one WD to another occurs on a dy- 
namical timescale if the mass ratio (donor/accretor) exceeds 
0.628 and the outcome is coalescence of the WDs. If two He 
WDs merge it is assumed that the temperature produced is 
enough to ignite the triple-Q reaction and that the nuclear 
energy released destroys the star. Four of the ten He-He WD 
binaries have mass ratios that satisfy this condition and two 
of them remain in the cluster at 4 Gyr. Iben (1990) discusses 
the possibility of subdwarf-0 or B stars forming from merged 
He WDs which ignited helium at the base of the accretion 
layer. This alternative scenario of helium star formation had 
previously been mentioned by Webbink (1984). 

6.3 Luminosity Functions 

The luminosity function (LF) for the 616 single MS stars in 
Model 2 at 4 Gyr is shown in Figure ITCl This is compared 
to the LF of the initial model (12 000 single MS stars) and 
also the LF at 4 Gyr for the model evolved with the pop- 
ulation synthesis code. The latter distribution we call the 
non-dynamical LF and it contains 11 341 stars up to the MS 
turn-off at V = 13. So 659 stars have evolved off the MS 
by 4 Gyr. The primordial and non-dynamical MS LFs are 
identical for V > 16 and very similar up to V = 13 apart 
from some fluctuation owing to the evolution of near-turn- 
off stars. So the non-dynamical MS LF at 4 Gyr can be used 
to infer the IMF of the population. If the slope of the clus- 
ter MS LF matches that of the non-dynamical LF it would 
also retain an imprint of the IMF. However, we can see from 
Figure 1161 that the two are far from being a good match. 
The cluster LF has been significantly flattened by the pref- 
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erential evaporation of low-mass stars from the cluster and 
by merging within binaries creating new MS stars near the 
turn-off. 



In Figure nn we split the cluster into five radial regions 
to examine the radial dependence of the MS LF slope and in- 
vestigate whether there is any region which retains sufficient 
information about the IMF. There is a clear radial depen- 
dence with the central region containing many more mas- 
sive stars than low-mass stars and the slope becomes flatter 
as we move outwards through the cluster. The normalized 
non-dynamical MS LF is compared to the cluster LF in each 
region and we can see that it is only in the outer regions that 
the slopes converge but, even then, the cluster LF is deficient 
in very low-mass stars and has an over- abundance of stars 
near the turn-off. Terlevich (1987) pointed out that in non- 
isolated clusters, heating by the galactic tide has the effect of 
making the velocity distribution isotropic and as such orbits 
ejected from the centre will tend to avoid returning there. 
This effect is in addition to mass segregation and helps in 
understanding the presence of massive MS stars in the outer 
regions. The LF for the five regions combined is also shown 
in Figure lTTI and is essentially flat with a clear preference for 
massive MS stars compared to the non-dynamical (or field) 
distribution. Fan et al. (1996) show that the overall MS LF 
for M67 is quite flat and that the centre of the cluster is 
deficient in low-mass stars. They also found a radial depen- 
dence for the mass function slope. Bonatto & Bica (2003) 
confirm the tendency for massive stars in M67 to be more 
centrally concentrated and that the central regions of the 
cluster show a turnover in the cluster luminosity function at 
low masses. They find that the halo of the cluster is enriched 
in low-meiss stars and their MS LF for the outer regions of 
M67 gives the best match to the expected LF slope for field 
stars (Kroupa, Tout & Gilmore 1993). 



We have presented LFs of single MS stars as this is the 
ideal situation to deal with if one wants to make inferences 
about the IMF of the stellar population. Evolved stars such 
as giants complicate matters because they suffer mass loss 
and their luminosity evolution is not as well understood, and 
binary star evolution is obviously more uncertain than that 
of single stars. When dealing with observed data the goal 
is the same but the process is not as straightforward. Re- 
moving evolved stars from the LF by inspection of the CMD 
is not too difficult. Binaries with high mass-ratios can also 
be removed using the method of CMD inspection. However, 
low mass ratio binaries are problematic and will lead to con- 
tamination of the LF. In Figure lTTI (sixth panel) we compare 
the slope of the LF for single MS stars to that of the same 
LF but with binaries with mass ratios less than 0.5 included 
as well. There are 532 single MS stars and 380 low-q bi- 
naries and the distributions have been normalized so that 
the slopes may be compared. We see that the inclusion of 
the low-q binaries does not affect the LF slope except at 
the low-mass end where binaries are biased towards high- 
q. Therefore, our model indicates that contamination of the 
MS LF by binaries is not necessarily a problem when using 
the shape of the distribution to infer the slope of the IMF. 



7 DISCUSSION AND SUMMARY 

A conclusion to draw from inspection of our two binary rich 
open cluster models is that in an open cluster, at least, one 
requires a substantial population of seed binaries capable 
of making BSs in order to explain the observed numbers of 
BSs. For Model 1 we used the birth period distribution sug- 
gested by Kroupa (1995b) to generate the orbital parameters 
of the 9 000 primordial binaries. Binary population synthe- 
sis told us to expect these binaries to produce just one BS 
after 4 Gyr of evolution and this would be via Case C mass 
transfer in a wide binary. After performing the A'^-body sim- 
ulation we indeed found one BS in the model at 4 Gyr. For 
Model 2 we included 12 000 primordial binaries and chose 
the initial separation of each binary from a flat logarithmic 
distribution with an upper limit of 50 au. In this case we 
expected 25 BSs at 4 Gyr with 75 per cent being produced 
from Case A mass transfer in short-period systems. What we 
found after 4 Gyr of cluster evolution was 20 BSs with seven 
of these produced from non-perturbed primordial binaries. 
So the expected BSs were depleted while others were cre- 
ated. A major difference between the primordial binary se- 
tups of Models 1 and 2 is that the frequency of short-period 
binaries is much less in Model 1 and no BSs are expected 
to be produced from Case A mass-transfer in these binaries. 
The formation rate of non-primordial short-period binaries 
is very small in open cluster simulations and when it does 
occur it is generally as a result of a binary-binary encounter 
involving at least one existing short-period binary. Also, in 
the moderate density conditions of an open cluster we do 
not see BS formation from direct collisions of single stars 
while relatively wide primordial binaries may be destroyed 
in encounters with other stars before reaching a Case B or 
C mass-transfer stage. These factors explain the absence of 
BSs in Model 1 and tell us that the cluster environment on 
its own is not efficient at producing BSs. It is important to 
have a population of short-period binaries that either form 
BSs directly or become involved in dynamical encounters 
with other cluster stars and binaries and subsequently pro- 
duce BSs. Our results indicate that period distributions such 
as Kroupa (1995b) and Duquennoy & Mayor (1991) that 
predict minimal BS formation from the standard Case A 
mass-transfer scenario are not suitable as initial conditions 
for open cluster binaries. 

In Model 2 we observed formation pathways to explain 
the full variety of BSs found in M67, single BSs, super-BSs 
and binary BSs with a range of period-eccentricity combi- 
nations. At 4 Gyr we found that 50 per cent of the BSs are 
single and, of those in binaries, all but one has an eccentric 
orbit. This is in contrast to the binary population synthesis 
expectation for the primordial binaries of the model which 
gives a mix of 75 per cent single and only circular binaries. 
This tells us that dynamical encounters within the cluster 
environment play an important role in defining the nature 
of the BS population but perhaps not in boosting actual BS 
numbers. Less than half of the BSs were formed from pri- 
mordial binaries that did not have their evolution altered 
in any way by dynamical encounters and the majority of 
these involved Case A mass transfer followed by coalescence 
to form a single BS. Other cases involved perturbations to 
the orbits of primordial binaries that induced mass transfer 
or, in some cases, a delay of mass-transfer after a primordial 
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binary became involved in a short-lived exchange encounter 
and emerged intact but with its orbital parameters altered. 
Mass transfer also led to BS formation in binaries created 
from exchange interactions. The alternative to BS forma- 
tion via mass transfer is the collision of two MS stars at 
pericentre in a highly eccentric binary. This was observed 
to occur in primordial binaries with mildly eccentric orbits 
that had their eccentricity increased by perturbations re- 
ceived from nearby stars or binaries, often after becoming 
involved in a three- or four-body hierarchy. Exchange in- 
teraction also created highly eccentric binaries in which the 
component stars eventually collided to form a BS. Mapelli 
ct al. (2004) showed that the radial distribution of BSs ob- 
served in the globular cluster 47 Tucanae is best explained 
by formation from primordial binaries (mass-transfer sce- 
nario) in the outer regions of the cluster and a mixture of 
exchange-induced collisions and primordial binary evolution 
in the core. So even though the stellar environment is very 
different in the core of 47 Tucanae than in M67 it seems 
that a mixture of BS formation scenarios is favourable in 
both. Davies, Piotto & De Angeli (2004) demonstrated that 
exchange interactions can be detrimental to BS formation 
via primordial binaries in massive globular clusters. They 
consider that exchanges produce binaries with increasingly 
more-massive components and this could lead to binaries 
that would have created a BS to be observed in a globular 
cluster now being replaced by binaries that form BSs much 
earlier. Their experiments showed that an average number 
of encounters per binary of about 5 to 10 was beneficial for 
BS production and that anything in excess of this would se- 
riously hamper the BS formation rate. This result is tuned 
to the age and turn-off mass of globular clusters and is also 
dependent on the distribution of binary separations but wc 
note that in our open cluster simulations only a few per cent 
of the binaries were involved in multiple exchange encoun- 
ters. 

Because a significant fraction of the BSs in our model 
are infiuenced in some way by the formation of triple and 
quadruple subsystems it would seem prudent to consider in- 
cluding a primordial population of these. Especially as triple 
and higher-order systems are observed in young open clus- 
ters (Schertl et al. 2003) and the field (Tokovinin 1997). 
Sandquist (2004) mentions nine possible triple systems in 
M67. Also, van den Berg et al. (2001) identified the BS 
S1082 in M67 (classified as single by Milone & Latham 1992) 
with a spectroscopic close binary and based on X-ray emis- 
sion it is now thought to be a triple system containing an 
RS CVn star with the BS as the outer component (van den 
Berg et al. 2004). Our M67 model at 4Gyr contains 20 hi- 
erarchical triples and there are 20 to 30 stable triples in 
the Model 2 simulation at any time. None of the BSs in the 
model at 4 Gyr are in triple systems although one BS-binary 
is loosely bound to another binary with an orbital separa- 
tion of 5 X 10* au. There are instances of triples containing 
a BS at other times, primarily as a result of a BS forming 
in an eccentric binary collision within a four-body or higher 
order system and remaining bound to other members of the 
system. These are not generally long-lived. Short-lived triple 
systems arc also found to be responsible for inducing an ec- 
centricity in the orbits of previously circular BS binaries. 
This is the Kozai effect (Kozai 1962) produced by a cyclic 
relation between the inner eccentricity and the orbital in- 



clination and which is modelled in NB0DY4 (Aarseth 2003). 
One of the RS CVn binaries in Model 2 at 4 Gyr (#1568) 
is the inner binary of a triple system with a l.SM© MS star 
and an outer period of about lOOOd. The capability to in- 
clude primordial triple and quadruple systems has recently 
been added to NB0DY4 and we plan to utilize this in future 
simulations. 

Mathieu et al. (2003) report the finding of two spectro- 
scopic binaries in M67 with a position in the CMD about 
1 mag below the subgiant branch. They are high-probability 
proper-motion members and have also been confirmed as X- 
ray sources (van den Berg, Verbunt & Mathieu 1999). Based 
on their CMD position the primaries of those binaries are 
termed sub-subgiants and it has been postulated that they 
arc the products of stellar encounters on non-standard evo- 
lutionary tracks (Mathieu et al. 2003). One binary is circular 
with a period of 2.82 d and the other has an eccentricity of 
0.21 and a period of 18.4 d. The eccentric binary is in the 
core of M67 and shows high reddening which implies a sub- 
giant with extinction (Mathieu et al. 2003). The star lying 
below the subgiant branch in our M67 model was formed in 
a CE event that saw 0.29Mq of material ejected from the 
star when the cores merged. Extinction from this material 
was not considered when calculating the CMD position of 
the star but this possibility and its position below the sub- 
giant branch make it a candidate sub-subgiant explanation. 
However, to explain the binary nature the star would need to 
have been a member of a triple system in which the merged 
star remained bound to the third component. This is not an 
unrealistic scenario and the inclusion of primordial triples 
would make it more likely. 

The M67 model we have presented here is a marked im- 
provement on the model previously reported in Hurley et al. 
(2001). A major reason for this is in the construction of the 
models. The current model is the result of an iV-body sim- 
ulation that was evolved from zero-age whereas the Hurley 
et al. (2001) A^-body model started at 2.5 Gyr. Also Hurley 
et al. (2001) subjected their model to an unusually strong 
tidal field in an attempt to produce a cluster that contained 
the observed mass of M67 within a tidal radius of 10 pc. 
More recent observations suggest that the tidal radius of 
M67 is actually greater than this. To achieve our preferred 
M67 model (Model 2 at 4 Gyr) we did make an alteration 
to the tidal field. Model 2 placed the cluster on an orbit 
at 8 kpc from the Galactic centre as opposed to 8.5 kpc for 
the standard orbit used in Model 1. This is well within the 
bounds allowed by the actual orbit of M67 and an orbital 
speed of 220kms~^ was used for both Models 1 and 2. Hur- 
ley et al. (2001) used 350kms~^. So the modification of the 
tidal field for Model 2 is nowhere near as extreme as for the 
Hurley ct al. (2001) model. In fact, according to Baumgardt 
(1998) it is the tidal radius at perigalacticon that deter- 
mines the cluster dissolution timescale and considering that 
M67 has a perigalacticon of 6.8 kpc wo could have chosen an 
even stronger tidal field for Model 2, and thus started with 
a greater rmmbcr of stars and binaries. 

Results are also better for our new M67 model because 
half of the BSs are found in binaries compared to only one 
out of 22 found by Hurley ct al. (2001). The raw rmmbers 
of BSs are similar but the Hurley et al. (2001) model actu- 
ally did better at creating BSs by dynamical means - they 
expected only half as many BSs as their model produced 
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whereas our model actually produces slightly less than ex- 
pected. Comparison of the two models is not straightfor- 
ward because different distributions were used for the pri- 
mordial binary separations and other factors such as the 
tidal field used by Hurley ct al. (2001) significantly alter 
the evolution of the cluster. Also, the semi-direct nature of 
the Hurley et aJ. (2001) model, which started with 5 000 
single stars and 5 000 binaries at an age of 2.5 Gyr, com- 
plicates direct comparison. Our model with initially 12 000 
single stars and 12 000 binaries actually contained 4 322 sin- 
gle stars and 4 859 binaries after 2.5 Gyr of evolution so the 
numbers are comparable at this stage. However, in an at- 
tempt to compensate for skipping 2.5 Gyr of A'^-body evolu- 
tion Hurley et al. (2001) selected their primordial binaries 
using a mass function biased towards higher masses and then 
evolved these to an age of 2.5 Gyr with the BSE code. Com- 
parison of this population with our 4 859 binaries at 2.5 Gyr 
shows that the Hurley et al. (2001) model has a higher pro- 
portion of binaries with a combined mass in excess of the 
cluster turn-off mass at 4 Gyr. The semi-direct model also 
evolved at higher central density (pc > 1 000 stars pc~^ for 
the majority of the simulation) and with a smaller half-mass 
radius. This would explain the increased incidence of dynam- 
ical BS formation. 

It is reassuring that the structure of our M67 model 
provides a good match to the observed structure of M67. 
The half-mass radius, the tidal radius and the mass of the 
model and cluster arc all in agreement. The surface density 
profile of the model is also able to reproduce the observed 
M67 profile, except in the core where the model has too 
many stars. We have constructed a surface brightness pro- 
file for the model using the same software as used for treating 
observed data and found that it is very well fitted by an El- 
son, Fall & Freeman (1987) distribution. This distribution 
was determined from young and intermediate age clusters in 
the Large Magellanic Clouds which do not show the same 
degree of tidal truncation as the old globular clusters of our 
Galaxy to which the King (1962) models were fitted. Dif- 
ferences between the Elson, Fall & Freeman (1987) model 
and the King (1962) model are only apparent near the tidal 
radius of a cluster. The luminosity functions for MS stars 
in various spatial regions of our M67 model also give good 
agreement with the observed behaviour. Our M67 model has 
done a very good job of simulating the blue straggler and 
RS CVn stars observed in M67. So the model places strong 
constraints on the parameters of the primordial binary pop- 
ulation in this cluster. Also, by matching the properties of 
populations for which observations are likely to be complete 
(BSs and RS CVn stars), we can use the model to make pre- 
dictions about the completeness of other populations (e.g. 
white dwarfs and BY Draconis stars). 

The process of tailoring a simulation to the parameters 
of a particular cluster and taking the time to compare the 
real and model data in a consistent manner has certainly 
proved fruitful. An effort along these lines was made previ- 
ously by Portegies Zwart et al. (2001) who looked at young 
open clusters such as the Hyades using models of AT ~ 3 000, 
whereas we have focussed on old clusters where the interac- 
tion between dynamical, stellar and binary evolution is more 
pronounced. We look forward to taking this approach with 
other rich open clusters of various ages, such as NGC 2099 
(Kalirai et al. 2001) and NGC 188 (Stetson, McClure & Van- 



denBerg 2004), in order to further understanding the ini- 
tial conditions, dynamical history and stellar populations of 
these objects. It can then be taken into the globular cluster 
regime when the petaflops speed GRAPE-DR (Makino et 
al. 2003) eventually becomes available. We also urge those 
working on observations of star clusters to take full advan- 
tage of dynamical models when interpreting data. 
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Figure 3. The mass profile of Model 2 at 4 Gyr. The solid line 
rciHCScnts the total cluster mass and the dashcd-dot line corre- 
sponds to the luminous mass. Dashed vertical lines show the tidal 
radius (15.2 pc) and the half-mass radius for the luminous mass 
(2.7 pc). 
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Figure 1. Evolution of binary fraction with age for an AT-body 
cluster model starting with 18 000 stars and 2 000 binaries (Shara 
& Hurley 2002: solid line). Also shown is a model that started 
with 12 000 stars and 8 000 binaries (Hurley & Shara 2003: dashed 
line). 
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Figure 4. Mass profiles for Model 2 at 4 Gyr scaled by the 
total mass involved in constructing the profile, all stars (solid 
line) , luminous mass (dashed-dot line) and the blue straggler mass 
(dotted line). 
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Figure 2. The mass profile of Model 1 at 4 Gyr. The solid line 

represents the total cluster mass and the dashed-dot line is the lu- 
minous mass. Dashed vertical lines show the tidal radius (20.8 pc) 
and the half-mass radius for the luminous mass (4.2 pc). 
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Figure 5. Evolution of the number density of stars witliin the core (solid line) and the half-mass radius (dashed line) for Model 2. 
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Figure 6. Model 2 at 4Gyr shown in the observed plane (the yz- or transformed YZ-plane). 
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Figure 7. Surface density profile of M67 from 2MASS data pro- 
vided by Bonatto (private communication: open circles). Profiles 
from the model for all stars (x symbols) and for only luminous 
stars with mass greater than 0.8Mq (open star symbols) are also 
shown. The distributions have not been normalized, la error bars 
have been included for the observed profile but for the sake of clar- 
ity we have not included error bars for the model points although 
the errors will be of comparable magnitude. 
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Figure 8. Surface brightness profile for Model 2 at 4Gyr (solid 
data points) constructed using software provided by Mackey & 
Gilmore (2003). The projection is along the Y-axis. Error bars are 
calculated according to the method outlined by Djorgovski (1988, 
in lAU Symp. 126, p333) where a given annulus is split into eight 
sectors of equal area. The internal error for the annulus is the 
standard deviation of the surface brightness values for the eight 
sectors. A fit to the data of an Elson, Fall & Freeman (1987) model 
(solid line) gives a core radius of 0.66 pc. This is a three-parameter 
model V{r) = Vb(l + /a^)~~'''^ where Vb is the central surface 
brightness and a is related to the core radius by rc = a{2'^^~' — 
1)^/-^. The error of the fit is 1.05 which represents the sum of 
the squares of the differences between the data and model points 
in each bin weighted by the error for that bin. 
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Figure 9. Same as Figure ISl but restricted to stars in the range 
12 < y < 17. The main-sequence turn-off is at V ~ 13 (cor- 
responding to a mass of 1.32Mq) so the range covers one mag- 
nitude above the turn-off, approximately half-way up the giant 
branch, and four magnitudes below the turn-off, down to a mass 
of 0.75Mq. The surface density profile fit has a scaled error of 
0.25 and gives a core radius of 0.64 pc. 
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Figure 10. Cluster colour-magnitude diagram at 4.0 Gyr. Note that all binaries are assumed to be unresolved. To convert the luminosities 

and effective temperatures to magnitude and colour we have used the bolomctric corrections given by Kurucz (1992) and, in the case of 
WDs, Bergeron, Wesemael & Beauchamp (1995). A distance modulus of 9.7 (Hurley et al. 2001) has been assumed to place the simulated 
cluster at the distance of M67. 
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Figure 11. As for Figure [TT7I but for stars within the half-mass radius (3.8 pc) of the cluster (left panel) and exterior to the half-mass 
radius (right panel). 
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Figure 12. Distribution of periods and eccentricities for binaries 
found to contain a BS that were present in tlic M67 simulation 
(Model 2) at an age of 2 Gyr or later. The solid stars represent the 
orbital periods at the time of formation, when one of the stars 
became a BS or when a BS was exchanged in to a new binary 
(31 points). If the orbital parameters of any of these binaries 
subsequently experienced a significant change (Ae > 0.05 and/or 
AlogP/d > 0.1), owing to binary evolution or a perturbation, 
this is represented by an open star (21 points). Orbital parameters 
for the six BS binaries known to reside in M67 are denoted by 
open circles. 
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Figure 13. Number of BSs (solid line), BS-binaries (dashed line) 
and circular BS-binaries with P < 100 d (dotted line) during the 
M67 simulation (Model 2). 
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Figure 14. X-ray luminosity as a function of orbital period for 
MS-MS binaries with a primary mass of I.OMq or less (open cir- 
cles) and for MS- WD binaries (solid triangles). The upper solid 
line shows the X-ray luminosity for a I.OMq MS star at an age of 
4 Gyr calculated from Equation 2 as a function of orbital period, 
assuming that the star is in a binary of that period and experi- 
encing synchronous rotation. The lower solid line is for a O.IM0 
MS star. We also include data points (solid stars) for the possible 
BY Diaconis systems identified by van den Berg et al. (2004) for 
wliicli an orbital period is known {Lx values for these systems 
supplied by M. van den Berg, private communication). 
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Figure 15. Mass fraction of white dwarfs in Model 2 at 4 Gyr as 
a function of radius. The radial bins arc chosen so that 150 stars 
are sampled in each bin. Results for all WDs (solid line) and for 
only single WDs and double WDs (dashed line) are shown. The 
dotted line at /wD = 0.1 is the value expected from isolated 
population synthesis of the initial stars. 
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Figure 16. Luminosity function of single main-sequence stars in 
Model 2 at time zero (dashed line, 12 000 stars) and at 4 Gyr (solid 
line, 616 stars). Also shown is the luminosity function for the 
initial stars evolved to 4 Gyr with the population synthesis code 
(dotted line, 11341 stars). The histograms are not normalized. 
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Figure 17. Luminosity function of single main-sequence stars for 13 < V < 23 in Model 2 at 4Gyr split into five radial regions. The 

regions arc chosen so that there is a similar immbcr of stars in each (about 106). The sixth panel (lower-right) shows the combined 
distribution. We compare this with the luminosity function for single MS stars and MS-MS binaries with q < 0.5 (dashed line). Also 
shown (dotted line with solid triangles) in each panel is the population synthesis MS luminosity function at 4 Gyr but with the number 
of stars in each bin reduced by a factor of 60 in order to aid comparison of the distribution slopes. 
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Table 1. Parameters of the starting models at time zero for simulations performed in this work 
(see text for details). 





Model 1 


Model 2 




9000 


12000 


Afb 


9000 


12000 




8.5 


8.0 


VQ/km/s 


220 


220 


density profile 


Plummer 


Plummer 


binary periods 


Kroupa 


flatlog 


Mo/Mq 


14405 


18687 


T,h,o/Myr 


300 


290 


Rt/pc 


34.4 


31.8 




4.3 


3.9 



Table 2. General results at 4Gyr for simulations performed in this work. 





Model 1 


Model 2 


M/Mq 


3175 


2037 


/b 


0.53 


0.60 




9 


13 


-Rt/pc 


20.8 


15.2 


Rh/pc 


4.9 


3.8 


Ml/M0 


1987 


1488 


Mlio/M© 


1730 


1342 


-Rh,Llo/PC 


3.0 


2.7 



Table 3. Stellar population results at 4 Gyr for simulations performed in this work (see text for 
details). 



Model 1 Model 2 





1 


20 


-'VBS.bin 


1 


9 


-'VBs/-'VMS,2to 


0.01 


0.18 


-Rh,Bs/pc 




1.1 




2 


6 


A^CV 


3 


1 


/wD 


0.16 


0.15 


Rh,WT>/P'^ 


0.3 


0.6 
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Table 4. Details of the blue stragglers observed in Model 2 at 4 Gyr. Columns are the star ID 
number for the BS, mass of the BS, its V-band magnitude and {B — V) colour (that of the BS or 
unresolved BS-binary), the radial position in the cluster at 4 Gyr, and the time at which the BS 
obtained its current mass. If the BS is in a binary the companion type is given in Column 7 followed 

by the companion mass, orbital period and eccentricity. The final column gives a classification of 
the evolution history of the BS using the following key: prim = primordial binary; A = Case A 
mass transfer leading to coalescence; B = Case B mass transfer; C = Case C mass transfer; coll 
= collision in eccentric binary; exch = exchange interaction; and pert = perturbation to orbit. 



ID# 


M/Mq 


V 


{B-V) 


r/pc 


To/Myr 


type 




P/d 


e 


history 


3289 


2.10 


10.95 


0.12 


1.03 


3613 


MS 


1.3 


81.3 


0.86 


exch— coll— coll 


1418 


2.09 


10.82 


0.20 


0.37 


3844 


giant 


0.7 


1.95 


0.00 


C-exch-B 


2203 


2.08 


10.94 


0.17 


1.10 


3480 


MS 


0.8 


363 


0.30 


pert-coll-exch-pert 


2411 


1.97 


11.61 


0.07 


0.49 


3657 










prim-pert-A 


2565 


1.89 


11.44 


0.22 


19.38 


3652 










prim-pcrt-coU 


1613 


1.88 


11.36 


0.34 


0.92 


2871 


MS 


0.9 


11749 


0.27 


cxch-A-exch 


2321 


1.88 


11.71 


0.13 


1.74 


3971 


MS 


0.7 


26.9 


0.65 


prim-pert— coll 


2737 


1.80 


11.56 


0.39 


0.91 


2549 










prim-A 


2855 


1.74 


11.84 


0.28 


2.11 


2946 










prim-A 


3835 


1.73 


11.90 


0.26 


0.31 


3798 


MS 


1.0 


19498 


0.69 


pert-coll-exch-coU 


2973 


1.69 


11.93 


0.38 


1.10 


3115 










prim-pert-coll 


3021 


1.67 


12.17 


0.27 


2.96 


3313 










prim-pert-A 


3157 


1.63 


12.04 


0.56 


3.08 


1948 










prim— A 


3121 


1.64 


12.32 


0.27 


0.66 


3885 


MS 


0.3 


28184 


0.72 


p crt-coU— exch 


3207 


1.61 


12.27 


0.34 


2.44 


3803 










prim-A 


3445 


1.53 


12.36 


0.60 


4.42 


2768 


MS 


0.3 


8511 


0.58 


pert-coU— exch 


3523 


1.51 


12.77 


0.35 


0.75 


3896 










prim— A 


3877 


1.40 


12.82 


0.52 


1.31 


1425 










prim— A 


3885 


1.40 


12.80 


0.54 


1.65 


1241 










prim— A 


1378 


1.36 


12.92 


0.57 


1.55 


1957 


WD 


0.6 


1660 


0.25 


prim-C-pert 
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Table 5. Detailed description of the formation scenario for some of the BSs listed in Table 



ID# explanation 

1378 This primordial binary began with an orbital period of 14454 d, an eccentricity of 0.83 and stellar masses of 1.75 and 1.21M0. 

After 1954 Myr the orbit began to circularize with the more massive star on the AGB. The primary filled its Roche lobe shortly 
afterwards when the orbit was circular with P = 2 203 d. At this point the masses were 1.02 and 1.23M0 so that CE evolution 
was avoided and stable Case C mass-transfer began. This phase ended when the primary had shed its entire envelope to become 
a O.62M0 CO WD in a circular binary of period 1 819 d with a 1.36M0 MS star companion. This was perturbed by a third star 
at T = 2 504 Myr when 0.35 pc from the cluster centre. The period was reduced to 1659d and an eccentricity of 0.25 was induced. 

1418 The proto-BS originated as a 1.23Af0 star in a binary of period 4 466d and eccentricity 0.7 with a 1.67Af0 companion. 

After 2 243 Myr the initially more massive star had evolved to the AGB and wind mass loss had reduced it to 1.21M0 while 
tidal forces had circularized the orbit (P = 1 750 d). At this point Case C mass-transfer began. This ended with a O.64M0 
CO WD and a 1.45M0 MS star in a circular orbit with P = 1445d. At T = 3 243 Myr it was involved in an exchange with 
a 1.35M0 star. The WD was ejected to leave a binary with e = 0.98 and P = 3631d. At T = 3 844 Myr the 1.35M0 star 
evolved off the MS, tides circularized the orbit and Case B mass-transfer started (ongoing at 4Gyr). 

1613 A primordial binary of period 47 860 d and eccentricity 0.8 containing stars of mass 1.46 and O.3M0 became involved 
in a four-body interaction with another primordial binary (P = 3 782d, e = 0.65 and masses of 0.42 and O.38M0) after 
2 089 Myr. From this a short-period eccentric binary containing the 1.46 and O.42M0 stars was formed (P = 1.1 d, 
e = 0.8). At T = 2 871 Myr the binary had circularized and the 1.46M0 MS star began Case A mass transfer quickly 
followed by coalescence. The I.88M0 single BS became involved in a three-body hierarchy at T = 2 932 Myr with a wide 
primordial binary (0.86 and O.8IM0) in the cluster core. Eventually the least massive star was ejected from the system. 

2203 The proto-BS began life as a 1.23M0 star with a O.85M0 companion in a primordial binary of period 363 d and an 
eccentricity of 0.3. After 3 479 Myr a binary-binary encounter left this binary as the inner component of a four-body 
system with stars of mass 0.83 and O.58M0. The eccentricity of the inner binary was driven up to 0.99 so that the 
stars collided and formed the 2.O8M0 BS. The O.58M0 star escaped the system and the BS remained bound to the 
O.83M0 MS star (P = 372 d, e = 0.84). A subsequent interaction reduced the period to 363 d and the eccentricity to 0.3. 

2321 This primordial binary was originally comprised of 1.29 and O.73M0 stars in a circular orbit with a period of 31 d. 

At T = 3 360 Myr while in the core of the cluster this binary had a close encounter with another primordial binary to 
form a quadruplet system. An eccentricity of 0.15 had been induced into the binary at this point. At T = 3 971 Myr 
the I.29M0 star collided with a O.59M0 MS star — the relative eccentricity of these two stars had reached 0.99 - to 
form the I.88M0 BS. The BS remained bound to its original companion and the fourth star (O.4IM0) was ejected. 

2411 This began in a circular primordial binary with Id period and component masses 1.29 and O.68Af0. Case A mass-transfer 
was expected to start at 4 320 Myr but perturbation at 3 100 Myr while binary was in the core hardened the orbit and 
mass-transfer began at 3 260 Myr. Angular momentum loss from the binary lead to coalescence of the stars at 3 657 Myr. 

2565 This started in a primordial binary with 3 236d period and eccentricity of 0.32. The orbit was wide enough that interaction 
between the 0.95 and O.94Af0 stars was not expected. Perturbation to the orbit while the binary was 1.0 pc from the cluster 
centre increased the eccentricity to 0.99. The orbit became chaotic and the stars collided and merged at T = 3 652 Myr. 

2973 This is similar to #2565. Interaction was not expected in a wide primordial binary with P = 4 075d and e = 0.59. A 

perturbation pumped the eccentricity to 0.99 so that the 0.89 and O.8OAf0 MS stars collided and merged at T = 3 115 Myr. 

3021 This originated in a circular primordial binary with P = 0.7d and component masses 1.01 and O.66M0 that was 

expected to begin Case A mass-transfer after 1 413 Myr. Prior to this, at T = 1 010 Myr, the binary was involved in a 
short-lived exchange encounter. The primordial binary emerged intact but the period had increased enough to delay 
the onset of mass-transfer until T = 2 900 Myr. The stars merged at T = 3 313 Myr. 

3121 This was a primordial binary composed of 1.09 and O.54M0 stars with a period of 562 d and an eccentricity of 0.31 which 
became part of a sextuplet after 3 762 Myr. The eccentricity of the binary was increased by perturbations from the 
other members until at T = 3 885 Myr the orbit became chaotic, the eccentricity reached 0.99 and the two stars collided. 
The resulting 1.64M0 BS remained bound to a 0.33Mq member of the sextuplet to give the binary observed at 4Gyr. 

3289 Two short-period circular primordial binaries (P = 1.3 d, Mi = O.82M0, M2 = 0.76Mq and P = 4.2 d. Mi = 1.25M0, 
M2 = O.5IM0) became embroiled in a four-body system after 3 611 Myr. The 0.82 and O.51M0 stars formed an inner 
binary and collided after the eccentricity was pumped up to unity. The collision product then collided with the O.76M0 
star to form a 2.O9Af0 BS. The BS remained bound to the Mi = 1.25M0 MS star in an eccentric orbit with P = 81d. 

3445 The proto-BS originated as a O.78M0 star in a 2.9 d circular orbit with a O.75M0 companion. No interaction between 
the stars was expected in this short-period primordial binary. After 2 768 Myr the binary formed a triple system with 
a O.25M0 star. The close presence of the third star induced an eccentricity of 0.84 into the binary orbit so that the 
stars came into contact and merged. The 1.53M0 merged star remained bound to the O.25M0 star. 

3835 A primordial binary consisting of 0.82 and O.6OM0 stars in a 1.8 d circular orbit became involved in a four-body system with 
another binary (formed earlier from an exchange) at T = 3 797 Myr. The stars in the primordial binary collided and merged 
and this new star then collided with a O.3IM0 star to form a 1.73A/0 BS. The BS remained bound to the fourth member of 
the system, a 1.O4M0 MS star, in a long-period eccentric orbit. 
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Table 6. Details of the RS CVn binaries observed in Model 2 at 4 Gyr. The first column gives the 
star ID number for the subgiant star and this is followed by the mass of the subgiant. Column 3 
gives the mass of the MS star companion. The V-band magnitude and {B — V) colour of the 
unresolved binary are then given. The sixth column gives the radial position in the cluster of the 

binary at 4 Gyr. The orbital period and eccentricity arc given in Columns 7 and 8, respectively. 
Some remarks on each binary are provided in the final column. 



ID# Mi/Mq M2/MQ V (B-V) r/pc P/d e remarks 



1568 


1.33 


0.95 


12.50 


0.60 


0.27 


1.6 


0.0 


exchange 


1799 


1.36 


1.08 


12.30 


0.88 


4.91 


5.0 


0.0 


primordial 


2335 


1.33 


0.68 


12.69 


0.68 


0.28 


3.0 


0.0 


primordial 


2383 


1.37 


0.62 


10.49 


1.17 


1.32 


19.8 


0.0 


primordial 


2633 


0.63 


1.19 


13.82 


0.64 


1.11 


0.4 


0.0 


primordial; semi-detached 


2873 


1.33 


0.40 


12.05 


0.63 


8.46 


6.2 


0.0 


primordial 



